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DOD  Prostate  Cancer  Research  Program 


Principal  Investigator:  Lei  Xina.  Ph.D. 


I.  INRTODCUCTION 

This  Idea  Award  (DAMD17-03-1-0023,  entitled  “Intensity  Modulated  Radiation  Treatment  of  Prostate 
Cancer  Guided  by  High  Field  MR  Spectroscopic  Imaging”)  was  awarded  to  the  principal  investigator  (PI)  for 
the  period  of  May  1,  2003 — April  30,  2006.  This  is  the  annual  report  for  the  first  funding  period  (May  1, 
2003  -April  30,  2004).  The  goal  of  this  project  is  to  establish  biologically  conformal  -  as  opposed  to 
anatomically  conformal  -  IMRT  as  a  viable  modality  through  integration  with  3T  MRSI  imaging  to  more 
effectively  kill  prostate  tumor  cells.  The  underlying  hypothesis  driving  this  work  is  that  the  MRSI-guided 
IMRT  will  provide  substantially  improved  dose  distributions  required  to  achieve  greater  local  tumor 
control  while  maintaining,  or  reducing,  complications  to  sensitive  structures.  The  specific  aims  of  the 
project  are:  (1)  To  establish  a  robust  procedure  for  registering  and  mapping  of  MR  spectroscopic  data  to 
CT/MRI  images  for  prostate  irradiation.  (2)  To  develop  an  inverse  planning  system  for  MRSI-guided 
IMRT  prostate  treatment  and  demonstrate  the  feasibility  of  concurrent  dose  escalation  to  intraprostatic 
lesion(s)  through  a  set  of  phantom  studies  and  at  least  two  previously  treated  prostate  cases  who  had 
undergone  CT/MRSI  scans.  Under  the  generous  support  from  the  U.S.  Army  Medical  Research  and 
Materiel  Command  (AMRMC)  ,  the  PI  has  contributed  significantly  to  prostate  cancer  research  by 
applying  physics  and  engineering  knowledge  to  prostate  cancer  research.  A  number  of  conference  abstracts 
and  refereed  papers  have  been  resulted  from  the  support.  The  preliminary  data  obtained  under  the  support 
of  the  grant  has  also  enabled  the  PI  to  start  new  research  initiatives  and  significantly  advanced  his 
academic  career.  In  this  report,  the  past  year’s  research  activities  if  the  PI  are  highlighted. 

II.RESEARCH  AND  ACCOMPLISHMENTS 

Prostate  Intensity  modulated  radiation  therapy  (IMRT)  is  an  image-guided  process  whose  success 
depends  critically  on  the  imaging  modality  used  for  treatment  planning  and  the  level  of  integration  of  the 
available  information.  In  current  clinical  practice,  IMRT  prostate  treatment  planning,  performed  under  the 
guidance  of  CT  or  MRI,  is  aimed  at  achieving  uniform  dose  distribution  to  the  entire  prostate  volume. 
While  tumor  biology  plays  a  crucial  role  in  the  treatment  outcome,  neither  MRI  nor  CT  provides  biological 
information^,  11,  13).  Recent  advancement  of  MR  spectroscopic  imaging  (MRSI)(2,  10,  14),  especially 
high  field  (3  Tesla)  MRSI,  now  makes  possible  to  effectively  distinguish  between  regions  of  cancer  and 
normal  prostatic  epithelium  with  clinically  acceptable  resolution  (~3mm).  Coupled  with  the  technical 
capability  of  IMRT  for  delivering  customized  3D  dose  distributions,  an  important  and  timely  question  to  ask 
is  whether  the  metabolic  information  derived  from  MRSI  can  be  used  for  more  objective  tumor  delineation 
and  for  guiding  IMRT  treatment  planning. 
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Toward  establishing  biologically  conformal  IMRT  for  prostate  cancer  treatment,  we  have  done 
substantial  work  in  the  past  year  in  integration  of  MRSI  imaging  into  IMRT  treatment  planning  and  in 
improving  the  existing  prostate  IMRT  planning  techniques.  The  research  is  sorted  into  four  different 
categories  and  summarized  below. 

Development  of  3T  MRSI  techniques :  We  have  successfully  implemented  single-  and  multi-voxel  and 
versions  of  J-resolved  MRS  sequences  in  order  to  better  detect  and  quantify  a  number  of  smaller  1H  MRS 
peaks,  including  those  from  glutamate  (Glu)  and  myo-Inositol  (ml).  The  detection  of  many  smaller 
metabolites  in  the  1H  spectrum  is  complicated  by  signal  overlap  due  to  both  a  small  dispersion  in  chemical 
shift  and  the  multiplet  structure  of  J-couple  resonances.  An  effective  remedy  is  to  use  two-dimensional  J- 
resolved  MRS  (2D-J-MRS)  where  chemical  shift  and  J-coupling  information  are  separated  into  two  different 
frequency  axis  known  as  fl  and  f2  respectively.  Closely  related  to  2D-J-MRS,  Constant-Time  PRESS  (CT- 
PRESS)  has  been  introduced  as  a  method  to  detect  coupled  resonances  with  effective  homonuclear 
decoupling  and  high  SNR .  A  PRESS  module  for  spatial  localization  is  followed  by  an  additional  refocusing 
pulse  whose  position  is  shifted  from  excitation  to  excitation  to  encode  chemical  shift  in  the  second  frequency 
dimension  (fl).  In  contrast  to  2D-J-MRS,  the  time  interval  between  excitation  and  data  acquisition  is  kept 
constant.  While  CT-PRESS  is  a  2D  experiment,  data  are  usually  display  as  a  ID  “diagonal”  spectrum 
obtained  by  integrating  the  2D  magnitude  spectrum  along  f2  within  a  limited  range  about  the  spectral 
diagonal.  As  we  have  demonstrated,  an  advantage  of  this  technique  is  that  the  evolution  time  can  be 
optimized  with  respect  to  the  coupling  constants  of  a  given  spin  system,  hence  increasing  the  SNR  for  a 
particular  metabolite.  In  vivo  data  demonstrated  the  robust  detection  of  Glu,  NAA,  Cho,  Cre,  and  ml 
resonances.  We  have  further  extended  our  CT-PRESS  acquisition  with  the  addition  of  spiral  readout 
gradients  for  spectroscopic  imaging.  In  vivo  results  demonstrate  single-slice  variable  density  spiral  MRSI 
CT-PRESS  with  excellent  spectral  quality  in  good  agreement  with  a  corresponding  single-voxel  acquisition. 
A  manuscript  on  in  vivo  detection  of  citrate  for  prostate  cancer  at  3T  has  been  submitted  to  Magnetic 
Resonance  Imaging  in  Medicine. 

Using  Biomechanical  Model  to  Register  Endorectal  Coil-Based  MRI/MRSI  with  Treatment  Planning 
CT  Images:  As  it  is  practiced  now,  radiation  therapy  treatment  planning  is  mostly  done  based  on  CT 
image  and  it  is  thus  required  to  map  any  new  imaging  information  onto  treatment  planning  CT  images. 
The  introduction  of  endorectal  surface  coils  significantly  improve  spatial  resolution  and  signal-to  noise 
ratio  (SNR)  of  prostate  MR  imaging  and  allows  evaluation  of  the  tumor  location,  tumor  volume,  capsular 
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penetration,  invasion  of  neurovascular  bundle,  and  seminal  vesicle  involvement,  which  is  crucial  for 
accurate  treatment  planning.  Endorectal-coil  based  MRSI  has  also  been  shown  effective  in  distinguishing 
between  areas  of  cancer  and  normal  prostatic  epithelium  through  differences  in  [choline+  creatinej/citrate 
ratio  (2,  5,  6,  14).  However,  the  use  of  endorectal  probe  inevitably  distorts  the  prostate  and  other  soft 
tissue  organs,  making  it  impossible  to  directly  fuse  the  acquired  image  data  onto  treatment  planning  CT. 
In  figure  1  we  show  the  difference  between  endorectal  coil-based  MRI  defined  and  CT-defined  prostate 
volume  (8).  In  order  to  fuse  MRI/MRSI  with  treatment  planning  CT,  it  is  needed  to  develop  an  effective 
deformable  image  registration  procedure.  Otherwise,  the  gain  from  the  use  of  the  state-or-the-art  imaging 
techniques  may  be  lost  due  to  the  inferior  performance  of  image  registration. 


Figure  1  Difference  between  endorectal  coil-based  MRI  defined  and  treatment  planning  CT-defined 
prostate  volumes. 

Image  registration  is  the  determination  of  a  geometrical  transformation  that  aligns  points  in  one 
view  of  an  object  with  corresponding  points  in  another  view  of  that  object  or  another  object.  Different 
methods  have  been  proposed  to  deal  with  the  registration  of  deformed  organs  but  the  optimal  approach  is 
still  illusive  in  many  cases,  especially  when  the  distorations  are  severe.  We  have  recently  applied  the  thin- 
plate  splines  (TPS)  method  (7)  to  deal  with  the  registration  of  endorectal  MRI/MRSI  and  verified  the 
accuracy  of  the  registration  by  phantom  and  patient  studies  (the  work  has  been  conditionally  accepted  for 
publication  in  Medical  Physics ).  For  this  purpose,  a  thin  plate  spline  (TPS)  transformation  first  introduced 
by  Bookstein  (1)  was  implemented  to  establish  voxel-to-voxel  correspondence  between  a  reference  image 
and  a floating  image  with  deformation.  The  idea  is  to  find  a  continuous  transformation  within  a  suitable 
Hilbert  space  to  minimize  the  landmark  difference  in  two  images.  The  detailed  description  of  the  TPS 
transformation  can  be  found  in  Bookstein’s  original  paper  (1).  To  access  the  quality  of  the  registration,  an 
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elastic  phantom  with  a  number  of  implanted  fiducial  markers  was  designed.  Radiographic  images  of  the 
phantom  were  obtained  before  and  after  a  series  of  intentionally  introduced  distortions.  After  mapping  the 
distorted  phantom  to  the  original  one,  the  displacements  of  the  implanted  markers  were  measured  with 
respect  to  their  ideal  positions  and  the  mean  error  was  calculated.  Phantom  studies  showed  that  using  the 
deformable  registration  method  the  mean  landmark  displacement  error  was  0.62  ±  0.39  mm  when  the 
distortion  was  of  the  order  of  23.07  mm.  A  deformable  model  seems  to  be  necessary  to  faithfully  map  the 
metabolic  information  onto  the  treatment  CT  images.  When  a  non-deformable  method  based  on  a  rigid- 
body  transformation  and  scaling  was  used  for  the  same  distortion,  the  mean  displacement  of  the  fiducials 
with  respect  to  their  actual  positions  was  found  to  be  as  large  as  12.95  ±  0.57  mm.  In  patient  studies,  CT 
images  of  two  prostate  patients  were  acquired,  followed  by  3  Tesla  (3T)  MR  images  with  a  rigid 
endorectal  coil.  For  both  patient  studies,  significantly  improved  registration  accuracy  was  achieved.  The 
prostate  centroid  position  displacement  was  0.58  ±  0.10  mm  and  the  coincidence  index  was  92.6  +  5.1  % 
when  a  TPS  transformation  was  used.  Different  from  the  non-deformable  approach,  the  TPS-based 
registration  accommodates  the  organ  distortion  and  enables  us  to  achieve  significantly  higher  MR/MRSI 
and  CT  image  registration  accuracy.  A  more  advanced  and  automated  finite  element  method  is  being 
developed  to  attack  the  problem  (12),  which  is  briefly  summarized  here. 

The  finite-element  calculation  was  performed  using  an  open-source  software  toolkit  named  the 
Insight  Toolkit  (ITK)  (Insight  Software  Consortium).  For  the  problem  of  multi-modality  deformable 
registration,  the  ITK  provides  an  implementation  of  finite-element  model  (FEM)-based  algorithm.  In  this 
implementation,  each  image  is  viewed  as  a  set  of  iso-intensity  contours.  The  main  idea  is  that  a  regular 
grid  of  forces  deform  an  image  by  pushing  the  contours  in  the  normal  direction.  The  orientation  and 
magnitude  of  the  displacement  is  derived  from  the  commonly  used  instantaneous  optical  flow  equation. 
The  FEM-based  deformable  registration  models  the  image  as  a  block  of  elastic  material  on  which  the 
forces  are  applied.  The  forces  are  computed  based  on  the  local  derivatives  of  image  metrics  on  points  of  a 
coarse  grid.  In  general,  subject  to  the  same  forces,  higher  elasticity  results  in  smaller  displacements,  and 
higher  density  results  also  in  smaller  displacements.  The  algorithm  is  initialized  with  a  deformation 
vectors  assessed  by  a  spline-based  deformable  registration  applied  on  a  mesh  of  automatically  generated 
points.  The  elasticity  and  tissue  density  used  in  our  calculation  are  from  the  literature  and  fine-tuning  of 
the  parameters  is  sometime  required  to  better  map  the  MRI  and  CT  images.  The  accuracy  of  a  FEM 
deformable  image  registration  method  is  evaluated  by  using  elastic  phantom  with  known  physical 
properties.  The  goodness  of  multi-modality  deformable  registration  is  quantitatively  assessed  using  a 
mutual  information  metric.  The  metric  is  voxel-based  and  no  control  point  or  manual  intervention  from 
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the  user  is  needed.  In  addition  to  the  mutual  information  assessment,  other  software  tools  are  also 
implemented  to  visually  assess  the  registration  results  and  to  facilitate  the  interaction  with  radiation 
oncologists.  Two  patient  studies  yield  consistent  agreement  and  suggest  that  the  approach  is  adequate  for 
radiation  therapy  application.  The  whole  registration  procedure  for  a  complete  3D  study  containing 
512x512x64  voxels  is  less  than  1  hour  on  a  standard  PC.  The  restored  MRI/MRSI  image  data  are  written 
in  the  CT-coordinates  in  DICOM  image  form  and  can  be  easily  imported  by  any  treatment  planning 
system.  The  registration  works  directly  on  3D  images  and  no  control  points  need  to  be  specified  by 
clinician.  We  are  in  the  process  of  refining  the  tool  for  MRSI-guided  prostate  IMRT  and  writing  up  the 
manuscript  for  publication. 


Figure  2.  A  photo  of  quality  assurance  phantom  built  for  testing/validating  the  geometric  and  metabolic  accuracy  of 
the  endorectal  coil-based  MRI/MRSI. 


Phantom  design  and  measurement  for  validation  and  quality  assurance  of  MRSI:  Increasingly,  ratios  of 
organ-specific  metabolite  intensities,  as  measured  by  MRSI  methods,  are  being  used  to  differentiate  regional 
levels  of  abnormality  for  diseased  tissue.  We  developed  a  QA  phantom  (figure  2  shows  a  photo  of  the 
phantom)  and  procedure  to  ensure  the  accuracy  of  MRSI-derived  metabolite  data  (3).  A  phantom  filled  with 
a  prostate-mimicking  solution  was  custom-built  with  an  insert  holding  a  few  vials  containing  calibrated 
solutions  of  precisely  varying  metabolite  concentrations  that  emulated  increasing  grade/density  of  prostate 
tumor.  This  work  is  an  extension  of  our  previous  work  on  brain  MRSI  phantom  measurement^).  Metabolite 
ratios  calculated  from  MRS  data  for  each  vial  were  compared  to  calibration  ratios  acquired  in  vitro  at  9.4  T. 
Least  squares  regression  analysis  was  used  to  investigate  these  relationships.  Regression  analysis  revealed  an 
expected  linear  relationship  with  regression  coefficients  close  to  1  and  intercepts  close  to  0  for  the  three 
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acquisition  modes.  This  simple  phantom-based  QA  of  MRS  metabolite  ratio  data  quality  is  useful  in  a 
routine  clinical  environment  or  for  development  of  new  MRS  and  MRSI  acquisition  software. 

Improve  IMRT  dose  distribution  by  using  non-uniform  penalty  scheme:  We  have  established  a  novel 
inverse  planning  formalism  for  producing  spatially  non-uniform  IMRT  dose  distributions.  A  strategy  to 
model  the  intra-structural  tradeoff  has  been  developed.  Generally  speaking,  voxels  within  a  structure 
volume  are  not  equivalent  in  achieving  their  dosimetric  goals.  To  consider  this  in  inverse  planning,  we 
proposed  to  purposely  modulating  the  penalty  on  individual  voxel  level  based  on  the  a  priori  dosimetric 
capability  (which  is  quantified  by  the  dosimetric  capability  that  measures  the  degree  for  a  voxel  to  achieve 
its  dosimetric  goal  and  is  obtained  by  using  a  Cimmino  algorithm  on  a  case  specific  basis).  We  have 
shown  that  an  inverse  planning  framework  with  customized  non-uniform  importance  factors  can  lead  to 
significantly  improved  IMRT  treatment  plans  that  would  otherwise  be  unattainable.  A  manuscript  has 
been  submitted  to  International  Journal  of  Radiation  Oncology,  Biology,  Physics  for  consideration  of 
publication. 

In  a  separate  but  related  project,  we  have  developed  a  biologically  more  sensible  yet  clinically 
practical  inverse  planning  system  for  functional  image-guided  IMRT  planning.  Current  inverse  planning  is 
done  using  dose-based  algorithms,  which  do  not  consider  the  nonlinear  dose  response  for  tumors  and 
normal  structures.  The  choice  of  structure  specific  importance  factors  represents  an  additional  degree  of 
freedom  of  the  system  and  makes  rigorous  optimization  intractable.  We  proposed  to  characterize  the  dose- 
volume  status  of  a  structure  by  using  the  concept  of  effective  volume  in  the  voxel  domain  and  constructed 
an  objective  function  with  incorporation  of  the  volumetric  information.  To  automate  the  determination  of 
the  structure  specific  importance  factors,  we  wrote  the  conventional  importance  factor  of  an  organ  into  a 
product  of  two  components:  (i)  a  generic  importance  that  parameterizes  the  relative  importance  of  the 
organs  in  the  ideal  situation  when  the  goals  for  all  the  organs  are  met;  (ii)  a  dose-dependent  factor  that 
quantifies  our  level  of  clinical/dosimetric  satisfaction  for  a  given  plan.  The  generic  importance  is 
determined  a  priori,  and  in  most  circumstances,  does  not  need  adjustment,  whereas  the  second  one,  which  is 
responsible  for  the  intractable  behavior  of  the  tradeoff  seen  in  conventional  inverse  planning,  is  determined 
automatically.  Compared  with  the  conventional  inverse  planning  technique,  we  found  that,  for  the  same 
target  dose  coverage,  the  critical  structure  sparing  was  substantially  improved  for  all  testing  cases. 

To  show  potential  benefit  of  the  stated  non-uniform  penalty  scheme,  in  Fig.  3  we  plotted  the  isodose 
distribution  for  a  prostate  IMRT  case.  With  this  newly  proposed  technique,  we  found  that  the  dose 
distribution  was  remarkably  improved  in  comparison  with  the  conventional  IMRT  plan  obtained  with 


9 


DOD  Prostate  Cancer  Research  Program  Principal  Investigator:  Lei  Xina.  Ph.D. 

structurally  uniform  importance  factors.  In  this  case,  the  dose  everywhere  within  the  sensitive  structure  was 
reduced  by  -17%.  Clinically,  this  implies  that  the  target  dose  may  be  escalated  by  -10%  while  keeping  the 
toxicity  at  the  level  achievable  by  current  IMRT  planning  technique.  Considering  the  improvement  was 
accomplished  purely  from  a  more  physical  modeling  of  the  system,  the  result  is  quite  striking.  The  above 
planning  formalism  is  ideally  suitable  for  deriving  spatially  non-uniform  dose  distributions.  This  work  is 
still  in  progress. 


Fig.  3  Comparison  of 
isodose  distributions  for  a 
5-field  prostate  case:  (a) 
conventional  approach 
(upper);  (b)  new  approach 
(lower).  The  contours  of 
the  target  (white)  and  the 
rectum  (green)  are  also 
shown.  It  is  clearly  seen 
that  the  use  of  non-uniform 
importance  factors  allow 
us  to  “push”  the  high  and 
intermediate  isodose 
curves  toward  the  target. 
The  dosimetric 

improvement  can  also  be 
seen  in  the  DVHs  shown 
in  Fig.  4. 


Fig.  4  Comparison  of  DVHs 
for  a  5-field  prostate  case 
(isodsoe  distributions  are 
shown  in  Fig.  3).  The  dashed 
(red)  and  solid  (blue)  curves 
correspond  to  conventional 
plan  and  the  plan  obtained 
using  non-uniform 

importance  factors, 

respectively. 


The  Idea  Award  for  Prostate  Cancer  Research  from  US  Amy  Medical  Research  and  Materiel 
Command  also  provides  a  unique  educational  opportunity  for  training  junior  researcher  through  the 
participation  of  research  activities.  In  this  aspect,  the  postdoctoral  fellow,  Eduard  Schreibmann,  has  been 
benefited  greatly  from  the  support.  After  obtaining  his  Ph.D.  degree  in  Medical  Physics  from  University  of 
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Patras,  Greece,  Dr.  Schreibmann  jointed  my  group  last  year.  Since  he  came  here,  he  has  had  opportunity  to 
learn  the  clinical  prostate  treatment  planning,  simulation,  and  quality  assurance  and  many  aspects  of 
radiation  treatment  of  prostate  cancer.  He  has  developed  an  elegant  and  clinically  practical  method  of  beam 
orientation  optimization  for  prostate  IMRT.  He  has  also  done  substantial  work  on  MRSI-CT  image 
registration,  which  will  be  presented  in  the  upcoming  ASTRO  annual  meeting.  He  will  continue  working 
under  the  Pi’s  supervision  on  the  project.  Given  his  training  obtained  under  the  support  of  this  grant  and 
his  performance  in  the  past  year,  I  expect  that  he  will  become  a  leading  researcher  in  the  fields  of  medical 
physics  and  prostate  radiotherapy  in  the  years  to  come. 


III.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Implemented  3T  prostate  MRSI  technique  with  J-resolved  pulse  sequence. 

•  Developed  a  TPS-based  deformable  image  registration  tool  for  mapping  prostate  MRSI  data  onto 
treatment  planning  CT/MRI  images.  A  more  automated  finite-element  deformable  registration 
method  to  further  improve  the  TPS  technique  is  close  to  be  completed. 

•  Validated  the  deformable  image  registration  model  using  elastic  phantom. 

•  Designed  and  fabricated  a  phantom  and  carried  out  a  series  of  measurement  to  validate  the  prostate 
MRSI  data  and  to  ensure  the  quality  of  endorectal  coil-based  MRI/MRSI  images. 

•  Developed  a  novel  inverse  planning  algorithm  with  spatially  non-uniform  penalty  for  optimizing  the 
prostate  IMRT  plan. 
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IV.  REPORTABLE  OUTCOMES 

The  following  is  a  list  of  publications  resulted  from  the  grant  support.  Copies  of  the  publication 
materials  are  enclosed  with  this  report. 

Refereed  publications: 

1 .  Yong  Y,  Xing  L,  “Inverse  planning  with  adaptively  evolving  voxel-dependent  penalty  scheme”, 
Medical  Physics ,  accepted. 

2.  Lian  J,  Hunjan  S,  Daniel  B,  Lo  A,  Levin  J,  Cardenas  C,  Dumoulin  C,  Watkins  R,  Rohling  K, 

Giaquinto  R,  Boyer  A,  Xing,  L.  Mapping  of  the  Prostate  in  Endorectal  Coil-Based  MRI/MRSI  and 
CT:  A  Deformable  Registration  and  Validation  Study.  Medical  Physics  31:  conditionally  accepted., 
2004. 

3.  Schreibmann  E  and  Xing  L,  “Beam  orientation  class-solutions  for  IMRT  prostate  cancer  treatment”. 
Medical  Physics,  submitted. 

4.  Yang  Y  and  Xing  L,  “Clinical  knowledge-based  inverse  treatment  planning”.  International  Journal 
of  Radiation  Oncology  Biology  Physics,  submitted. 

5.  Kim  D,  Mayer  D.,  Xing  L,  Daniel  B,  Margolis,  D.,  Spielman  D.,  “In  vivo  detection  of  citrate  for 
prostate  cancer  at  3  Tesla”,  Magnetic  Resonance  Imaging  in  Medicine ,  submitted. 

Conference  abstracts: 

The  Pi's  group  has  also  been  active  in  disseminating  our  research  results.  The  following  are  some  of 
the  presentations  given  in  various  national/international  meetings. 

1.  Lian  J.  and  Xing  L.,  Biological  Model  Based  IMRT  Optimization  with  Inclusion  of  Parameter 
Uncertainty,  International  Conference  on  Computers  in  Radiotherapy,  Seoul,  Korea,  May,  2004. 

2.  Shou  Z.  and  Xing  L.,  Improve  IMRT  dose  distribution  using  spatially  non-uniform  importance  factors. 
International  Conference  on  Computers  in  Radiotherapy,  Seoul,  Korea,  May,  2004. 

3.  Xing  L  and  Shou  Z.,  “Intrinsic  spatial  heterogeneity  in  inverse  planning  and  its  significant  role  in 
defining  the  IMRT  solution  space”,  2004  ASTRO  Annual  Meeting. 

4.  Yong  Y.  and  Xing  L,  “Clinical  knowledge  -based  inverse  planning”,  2004  ASTRO  Annual  Meeting. 

5.  Schreibmann  E.  and  Xing  L.,  “Dose-Volume  Based  Ranking  of  Incident  Beams  and  its  Utility  in 
Facilitating  IMRT  Beam  Placement”,  2004  ASTRO  Annual  Meeting. 

6.  Schreibmann  E.,  D.  Kim,  S.  L.Hancock,  A.L.  Boyer,  D.  Spielman,  B.  Daniel,  and  Xing  L.,  “Using 
Finite-Element  Method  to  Register  Endorectal  Coil-Based  MRI/MRSI  with  Treatment  Planning  CT 
Images”,  2004  ASTRO  Annual  Meeting. 

Invited  talks: 

The  PI  has  been  invited  as  abstract  review  and  session  chair  in  several  national/intemational  meetings 
and  has  given  a  number  of  invited  talks.  The  following  is  a  partial  list  of  the  relevant  activities  in  last  year. 

1 .  Debate  Session:  IMRT  or  Proton?  World  Congress  on  Medical  Physics  and  Biomedical  Engineering, 
Sydney,  Australia,  Aug.,  2003. 

2.  Dose  Matching  in  IMRT,  World  Congress  on  Medical  Physics  and  Biomedical  Engineering,  Sydney, 
Australia,  Aug.,  2003. 

3.  Integration  of  Functional  Imaging  into  Radiation  Therapy  Planning,  Workshop  on  Image  Guided 
Radiotherapy,  Lake  Tahoe,  Sept.,  2003. 

4.  Functional  and  Molecular  Imaging,  Midwinter  Symposium  of  AAPM  Southern  California  Chapter, 
Los  Angeles,  Jan,  2004. 

5.  Clinical  Knowledge-Based  Inverse  Planning,  International  Workshop  on  Intensity-Modulated 
Radiation  Therapy,  Medical  Imaging,  and  Optimization,  University  of  Haifa,  Israel,  June,  2004,  to  be 
delivered. 
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6.  Recent  Progress  in  Inverse  Treatment  Planning,  Annual  Radiological  Society  Meeting  of  Taiwan, 
Taipei,  March,  2004. 

7.  Functional  and  Molecular  Imaging,  Image-Guided  Interventions:  Trial  Design  Methodology 
Workshop,  Sponsored  by  the  Cancer  Imaging  Program,  NCI,  Feb.,  2004. 

8.  Toward  Biologically  Conformal  Radiation  Therapy.  International  Conference  on  Computers  in 
Radiotherapy,  May,  Seoul,  Korea,  2004. 

9.  Biological  Imaging  and  Biological  Imaging-Guided  IMRT,  American  Cancer  Society  El  Cerrito  Relay 
for  Life  Open  Ceremony,  May,  2004. 

Book  Chapters: 

1 .  Xing  L.,  Yang  Y.,  Spielman  D.,  Molecular/Functional  Image-Guided  Radiation  Therapy,  T.  Bortfetld, 
R.  Schmidt-Ullrich,  W.  de  Neve  (editors),  Spinger-Verlag  Heidelberg,  Berlin,  in  press. 


V.  CONCLUSIONS 

MR  Spectroscopic  Imaging-guided  IMRT  technique  is  being  developed  for  the  treatment  of  prostate 
cancer.  A  few  important  milestones  have  been  achieved  toward  the  general  goal  of  the  project.  These 
include  (i)  implemented  3T  prostate  MRSI  technique  with  J-resolved  pulse  sequence;  (ii)  development  of  a 
robust  image  registration  procedure  to  map  endorectal  coil-based  MRSI  data  onto  treatment  planning 
CT/MRI  images;  (iii)  experimental  validation  of  the  image  registration  tools;  (iv)  phantom  measurement  of 
MRSI;  (v)  improvement  of  currently  available  inverse  planning  algorithm  for  MRSI-guided  prostate 
IMRT.  Integration  and  further  refinement  of  the  above  tools  are  underway.  It  is  expected  these  tools  will 
greatly  facilitate  the  planning,  delivery,  and  quality  assurance  of  the  MRSI-guided  prostate  treatment. 
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Abstract 


In  current  inverse  planning  algorithms  it  is  common  to  treat  all  voxels  within  a 
target  or  sensitive  structure  equally  and  use  structure  specific  prescriptions  and  weighting 
factors  as  system  parameters.  In  reality,  the  voxels  within  a  structure  are  not  identical  in 
complying  with  their  dosimetric  goals  and  there  exists  strong  intra-structural  competition. 
Inverse  planning  objective  function  should  not  only  balance  the  competing  objectives  of 
different  structures  but  also  that  of  the  individual  voxels  in  various  structures.  In  this  work 
we  propose  to  model  the  intra-structural  tradeoff  through  the  modulation  of  voxel 
dependent  importance  factors  and  deal  with  the  challenging  problem  of  how  to  obtain  a 
sensible  set  of  importance  factors  with  a  manageable  amount  of  computing.  Instead  of 
letting  the  values  of  voxel  dependent  importance  to  vary  freely  during  the  search  process, 
an  adaptive  algorithm,  in  which  the  importance  factors  were  tied  to  the  local  radiation 
doses  through  a  heuristically  constructed  relation,  was  developed.  The  new  planning  tool 
was  applied  to  study  a  hypothetical  phantom  case  and  a  prostate  case.  Comparison  of  the 
results  w  ith  t  hat  o  btained  u  sing  c  onventional  i  nverse  p  lanning  t  echnique  w  ith  s  tructure 
specific  importance  factors  indicated  that  the  dose  distributions  from  the  conventional 
inverse  planning  are  at  best  sub-optimal  and  can  be  significantly  improved  with  the  help 
of  the  proposed  non-uniform  penalty  scheme. 


Key  word:  IMRT,  inverse  planning,  dose  optimization,  objective  function,  tradeoff. 
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Introduction 

Inverse  treatment  planning  is  an  indispensable  step  in  the  implementation  of  IMRT  and  its 
performance  critically  determines  the  success  of  an  IMRT  treatments'll.  While  the 
general  concept  of  inverse  planning  seems  to  be  logical  and  straightforward,  much 
remains  to  be  done  to  improve  the  currently  available  algorithms^*  12-14  Qne  of  the  key 
issues  that  have  been  overlooked  in  the  formulation  of  inverse  planning  problem  is  the 
intra-structural  tradeoff.  Indeed,  in  most,  if  not  all,  inverse  planning  algorithms,  the  dose 
prescription  and  weighting  factors  are  specified  on  a  structure  specific  basis.  An  implicit 
assumption  made  in  these  implementations  is  that  all  points  within  a  structure  are 
equivalent.  In  reality,  voxels  within  a  structure  are  generally  not  equivalent  in  complying 
with  their  dosimetric  requirements  and  this  inherent  heterogeneity  has  not  been  considered 
so  far.  Depending  on  the  patient’s  geometry,  beam  modality  and  field  configuration,  some 
regions  in  a  patient  may  have  better  chance  to  meet  the  prescription  than  others. 
Furthermore,  the  doses  at  different  regions  are  often  incompatible  and  even  compete  each 
other.  Appropriate  tradeoff  between  the  voxels  is  a  requisition  to  fully  exploit  the  potential 
of  IMRT  in  any  type  of  inverse  planning  techniques!  5-1 9 

In  this  work  we  propose  to  model  the  intra-structural  tradeoff  through  the 
modulation  of  voxel  dependent  importance  factors.  The  new  inverse  planning  technique 
allows  us  to  balance  not  only  the  competing  objectives  of  different  structures  but  also  that 
of  the  individual  voxels  within  any  structure.  Computationally,  the  determination  of  intra- 
structural  tradeoff  is  an  intensive  task  because  of  the  coupling  between  the  beamlet 
weights  and  local  importance.  To  obtain  an  adequate  set  of  local  importance  factors  with  a 
manageable  amount  of  computing  time,  we  develop  an  adaptive  algorithm.  The  main  idea 
of  the  adaptive  approach  is  to  heuristically  link  the  values  of  the  voxel-dependent 
importance  factors  with  the  corresponding  radiation  doses  and  to  continuously  update 
their  values  until  the  optimality  criterion  is  met.  The  new  planning  technique  has  been 
applied  to  study  a  hypothetical  phantom  case  and  a  prostate  case  and  its  advantage  over 
the  conventional  inverse  planning  technique  with  structurally  uniform  importance  has 
been  demonstrated.  The  algorithm  seems  to  afford  an  effective  way  of  modeling  the 
system  with  only  a  little  computational  overhead  and  shows  significant  promise  to 
improve  the  existing  inverse  planning  techniques. 
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Methods  and  Materials 


Theoretical  background 

Clinical  objectives  are  multifaceted  and  often  incompatible  with  one  another.  The  inverse 
planning  is  thus  inherently  a  multi-objective  problem:  each  anatomical  structure,  or  more 
generally,  each  voxel  within  a  structure,  as  an  independent  entity  with  its  own  treatment 
objective.  A  commonly  used  approach  is  to  combine  these  objectives  to  form  an  overall 
objective  function  by  weighted  sum  of  the  individual  objectives.  That  is 

f-E'A«A(o»,  (i) 

(7 —1 

where  the  summation  is  over  all  targets  and  sensitive  structures,  Dc(i)  is  the  calculated 
dose  in  voxel  i,  and  ra  and  Fa  are  the  weighting  factor  (or  the  importance  factor)  and  the 
objective  of  the  structure  indexed  by  o,  respectively.  Regardless  of  the  model  used, 
Fc  can  be  expressed  as  a  function  of  the  dose  distribution,  {Dc(i)},  within  the  structure. 
For  a  dose-based  optimization,  the  quadratic  form^O  21,  22 f 

F*=Trf,nmi)-D,(0]2,  P) 

™  cr  /=1 

is  often  used,  where  Na  represents  the  total  number  of  voxels  in  structure  a,  and  D0  (/)  is 

the  prescription  dose  in  voxel  i.  In  the  above  description,  two  types  of  weighting  factors 
are  involved:  the  inter-structural  importance,  {ra},  and  the  intra-structural  ones,  {r,}.  Up 
to  this  point,  however,  only  structure  specific  importance  factors  are  utilized  in  inverse 
planning  formulation  and  the  importance  of  voxel  specific  weighting  factors  has  been 
overlooked  other  than  the  fact  that  they  have  been  used  to  “tweak”  IMRT  dose 
distributions^"^, 18  indeed,  in  most,  if  not  all,  objective  functions,  the  voxels  within  a 
structure  are  tacitly  assumed  to  be  equivalent  in  complying  with  their  dosimetric 
requirements  and  their  relative  importance  factors  are  set  to  unity.  This  type  of  penalty 
scheme  seriously  limits  the  solution  space  and  often  leads  to  sub-optimal  plan.  To  give  a 
comprehensive  example,  one  can  imagine  the  consequence  when  two  or  more  structures 
in  a  system  are  restricted  to  take  a  single  importance.  To  be  able  to  assess  more  candidate 
plans  and  obtain  truly  optimal  IMRT  plans,  it  is  necessary  to  establish  a  voxel  dependent 
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penalty  scheme  in  which  the  penalty  at  a  voxel  depends  not  only  on  the  dose  deviation 
from  t  he  p  rescription  but  a  Iso  o  n  o  ther  p  hysical  a  nd  clinical  requirements  a  t  the  p  oint, 
which  can  be  reflected  by  purposely  modulating  the  {r,}. 

Voxel  dependent  importance  factor  and  its  heuristic  relation  to  the  local  dose 
While  the  introduction  of  voxel  dependent  importance  factors  has  the  advantage  for  us  to 
obtain  IMRT  solutions  that  would  be  otherwise  not  accessible,  a  practical  question  is  how 
to  determine  the  optimal  importance  distribution  fora  given  structure.  In  reality,  there 
may  be  more  than  one  source  that  contribute  to  the  voxel  heterogeneity  (for  example,  the 
in-equivalence  of  the  voxels  in  complying  their  dosimetric  goals  as  discussed  in  this  work, 
or  biological  heterogeneity  resulting  from  non-uniform  cell  density  or  radiation  sensitivity 
distribution)  and  the  local  importance  can  be  generally  written  as  a  product  of  the 
contributions  from  various  sources,  that  is,  rt  =  r'rfr,3...,  where  the  superscripts  index  the 

different  sources  of  contributions.  A  possible  approach  for  obtaining  the  importance 
distribution  is  so  called  a  priori  technique,  in  which  one  attempts  to  identify  the  origin  of 
voxel  in-equivalence  and  then  derives  the  spatial  distribution  of  {r,}  based  on  physical  or 
clinical  considerations.  One  can  also  proceed  in  a  posteriori  fashion  15-17  by  directly 
optimizing  the  {/*,},  similar  to  the  auto-selection  of  the  structure  specific  importance 
factors  proposed  by  Xing  et  all  7.  An  optimization  of  the  intra-structural  tradeoff  is  likely 
to  be  computationally  prohibitive  due  to  the  enormous  size  of  the  search  space  arising 
from  the  coupling  of  the  beamlet  weights  and  the  local  importance  factor.  To  circumvent 
the  problem,  in  this  work  we  develop  an  alternative  technique  in  which  the  importance 
distribution  is  determined  adaptively  under  the  guidance  of  a  heuristically  constructed 
function. 

Intuitively,  it  is  not  difficult  to  conceive  that  the  n  in  equation  (2)  can  be  regarded 
as  a  dose  dependent  parameter.  For  a  voxel  in  a  sensitive  structure,  for  instance,  the  value 
of  r,  should  be  higher  if  the  voxel  receives  a  high  dose  so  that  the  voxel  gets  more  penalty, 
and  vice  versa.  This  suggests  that  the  local  importance  can  be  expressed  as  a  function  of 
the  local  dose  and  updated  at  each  step  of  the  iterative  dose  optimization  process  without 
invoking  any  additional  computation.  While  the  general  monotonic  dependence  of  r,  on 
the  dose  is  clear,  its  specific  form  is  a  matter  of  experimentation.  We  propose 
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r,  =aa+alDc(i)k  (3) 

for  the  application,  where  a0,  ai,  and  kare  heuristic  structure  specific  parameters.  The 
value  of  k  is  greater  than  zero  for  a  voxel  in  a  sensitive  structure.  For  a  voxel  in  a  target,  k 
should  be  less  than  zero  in  order  to  assimilate  the  clinical  preference  over  overdosing 
versus  underdosing.  The  penalty  function  for  a  structure  consists  essentially  two  terms. 
The  first  term  is  the  conventional  one  with  uniform  importance  across  the  whole  structure 
and  the  second  term  modulates  the  importance  according  to  the  local  value  of  dose.  This 
phenomenological  importance  emphasizes  on  penalizing  those  voxels  receiving  a  high 
radiation  dose  yet  does  not  ignore  those  voxels  with  intermediate  and  even  low  dose  of 
radiation  at  each  step  of  the  iterative  calculation.  In  this  study,  we  set  cio  as  unity  and  a\  as 
where  Dref(i)  is  set  as  the  5%  tolerance  dose  TD5/5  for  sensitive  structures  and 

prescription  dose  D0  (/)  for  target,  k  is  an  adjustable  parameter  determined  by  a  trial-and 

error  process.  We  emphasize  that  the  {r,}  distribution  so  obtained  may  not  be  the  best 
possible  distribution  due  to  the  restriction  on  the  feasible  values  imposed  by  equation  (3). 
However,  as  will  be  seen  in  the  Results  section,  the  non-uniform  importance  factors 
resulted  from  the  use  of  equation  (3)  allow  us  to  significantly  improve  the  IMRT  dose 
distributions  obtained  using  the  existing  inverse  planning  algorithm  based  on  structure 
specific  importance  factors.  The  technique  thus  provides  a  practical  solution  to  what 
appears  to  be  a  computationally  overwhelming  problem. 

Calculation  process 

An  independent  optimization  module  based  on  the  objective  function  (2)  (with  {r,} 
determined  by  equation  (3))  was  integrated  into  the  PLUNC  treatment  planning  system 
(University  of  North  Carolina,  Chapel  Hill,  NC).  The  dose  calculation  engine  and 
varieties  of  evaluation  tools  existing  in  the  PLUNC  system  were  employed  for  this  study. 
For  a  given  patient,  a  set  of  structure  specific  importance  factors,  {ra},  were  determined 
empirically,  similar  to  that  in  the  conventional  inverse  planning.  For  a  given  set  of  {rff}, 

we  used  the  ray-by-ray  iterative  algorithm  (SIITP)  reported  in  References2^  24  t0  obtain 
the  optimal  beam  intensity  profiles.  The  only  difference  here  is  that  a  r,-modulated 
quadratic  objective  function  was  used.  The  calculation  time  for  a  typical  prostate  case  is 
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less  than  5  minutes  on  a  PC  with  P4  1.7GHz  and  1024MB  RAM. 

Case  Studies 

A  hypothetical  phantom  case  and  a  previously  treated  clinical  prostate  case  were  used  to 
evaluate  the  level  of  improvement  resulted  from  the  new  algorithm  with  adaptively 
determined  spatially  non-uniform  importance  factors.  The  results  for  both  cases  were 
compared  with  that  obtained  using  the  conventional  dose  optimization  method  with 
structural  specific  importance  factors.  For  the  purpose  of  plan  comparison  and  assessment, 
we  purposely  adjusted  the  structure  specific  importance  factors  in  such  a  way  that  the  dose 
coverage  of  the  target  volume  was  similar  for  the  two  IMRT  plans. 

The  phantom  case  consisted  of  a  target  and  a  sensitive  structure  (figure  la).  Five 
incident  beams  were  used  for  target  irradiation  (32°,  104°,  176°,  248°  and  320°  in  IEC 
convention).  The  incident  photon  energy  is  15  MY.  Table  I  listed  the  optimization 
parameters  for  both  the  approaches  for  this  case.  The  target  was  prescribed  to  100 
(arbitrary  unit)  and  the  prescription  dose  to  the  sensitive  structure,  D0(i)  in  equation  (2), 
was  set  to  zero. 

In  the  prostate  IMRT  case,  the  target  volume  was  the  prostate.  The  sensitive 
structures  involved  in  this  study  were  rectum,  and  bladder.  Five  equally  spaced  15MV 
photon  beams  (gantry  angles  of  0°,  72°,  144°,  216°,  and  288°  in  IEC  convention)  were 
used  for  the  treatment.  A  radiation  dose  of70Gy  was  prescribed  to  cover  99%  of  the 
target  volume.  For  the  sensitive  strictures,  the  prescription  doses  to  the  sensitive 
structures,  D0(i)  in  equation  (2),  were  set  to  zero.  The  optimization  parameters  for  both 
approaches  were  summarized  in  Table  II. 


Results  and  Discussion 

Hypothetical  phantom  IMRT  plans 

Figures  1  and  2  summarized  the  calculation  results  and  the  comparison  of  the  two  IMRT 
plans  obtained  using  the  newly  proposed  and  conventional  penalty  schemes.  In  figure  la 
we  depict  the  geometry  of  the  phantom  and  in  figures  lb  and  lc  we  show  the  isodose 
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distributions  for  the  plans  without  and  with  the  use  of  non-uniform  intra-structural 
importance  factors.  The  DVHs  of  the  target  and  sensitive  structure  are  plotted  in  figure  2. 
The  dashed  and  solid  curves  in  figure  2  represent  the  results  of  the  conventional  and  the 
newly  proposed  inverse  planning  approaches,  respectively.  It  is  clearly  seen  that  the 
proposed  technique  greatly  improves  the  critical  structure  sparing.  The  maximum  dose  to 
the  sensitive  structure  is  reduced  from  67  to  61  and  the  doses  to  other  voxels  in  the 
sensitive  structure  are  all  reduced  by  2-14%.  The  fractional  volume  reduction  in  the 
intermediate  and  low  doses  is  more  distinct.  For  example,  the  fraction  volume  receiving 
dose  above  40  is  decreased  from  19%  to  10%  while  that  receiving  dose  above  20  is 
decreased  from  95%  to  51%.  For  target,  the  volume  receiving  dose  above  85  is  similar 
for  the  two  algorithms,  98.4%  and  98.6%  for  the  conventional  and  the  new  approaches, 
respectively.  However,  the  volume  receiving  high  doses  is  decreased  (e.g.,  the  fractional 
volume  receiving  a  dose  above  102  is  reduced  from  20.7%  to  8.2%).  Interestingly,  the 
irradiation  to  the  normal  tissue  surrounding  the  target  and  sensitive  structure  seems  to  be 
less  in  the  new  optimization  technique.  For  example,  the  fraction  volumes  of  the  normal 
tissue  that  received  doses  above  60  and  40  are  decreased  from  12%  to  8%  and  from  34% 
to  25%,  respectively.  We  attribute  the  simultaneous  improvement  in  the  doses  of  all 
structures  to  the  greatly  enlarged  solution  space  when  non-uniform  importance  is 
permissible. 

Prostate  IMRT  plans 

The  IMRT  plans  obtained  using  the  two  different  approaches  for  the  prostate  case  are 
summarized  in  figures  3  and  4.  Figure  3  compares  the  isodose  distributions  in  two 
transverse  slices  and  one  sagittal  slice  for  the  two  plans.  The  DVHs  of  the  target  and  the 
sensitive  structures  for  the  two  plans  are  plotted  in  figure  4,  in  which  the  dashed  lines 
represent  the  results  obtained  using  the  conventional  algorithm  and  the  solid  lines  the 
newly  proposed  approach.  Similar  to  the  phantom  case,  when  non-uniform  importance 
factors  are  permissible,  the  conformality  of  the  final  dose  distribution  is  greatly  improved 
in  comparison  with  the  conventional  IMRT  planning  with  structurally  uniform 
importance  factors.  In  particular,  the  doses  to  the  sensitive  structures  are  dramatically 
improved.  As  canbe  seen  from  figure  3 ,  much  steeper  dose  gradients  atthe  interface 


8 


between  the  target  and  the  sensitive  structures  are  achieved.  From  the  DVHs,  it  is 
observed  that  the  rectum  and  bladder  are  better  spared.  For  example,  the  fractional 
volume  of  the  rectum  that  receives  a  dose  above  40Gy  is  dropped  from  36.5%  to  10.5%. 
Similarly,  10%  reduction  (from  30.5%  to  20.5%)  is  noticed  in  the  fractional  bladder 
volume  receiving  a  dose  above  30Gy.  The  integral  dose  to  the  normal  tissue  is,  once 
again,  reduced  as  a  consequence  of  the  use  of  non-uniform  importance  factors.  It  is 
important  to  emphasize  that  the  improvement  in  sensitive  structure  sparing  is  achieved 
with  essentially  no  change  in  the  dose  coverage  of  the  tumor  target.  For  the  target,  the 
fractional  volume  receiving  dose  below  70  Gy  is  decreased  from  about  2.1%  to  1.2%  and 
the  fractional  volume  receiving  dose  above  80  Gy  is  decreased  from  2.4%  to  1.5%.  The 
maximum  doses  for  the  two  plans  are  almost  same,  81.2Gy  and  81.4Gy,  respectively.  But 
the  minimum  dose  is  increased  from  63.4Gy  for  the  conventional  inverse  planning  to 
65.5Gy  for  the  new  technique. 

Finally,  it  is  interesting  to  point  out  that  the  EUD-based  optimization25-27  js  a 
special  case  of  the  empirical  model  proposed  in  this  work.  If  we  ignore  the  contribution 
from  the  dose  discrepancy  penalty  part  in  Eq.  (2)  (i.e.,  change  the  power  of  the  dose 
discrepancy  in  Eq.  (2)  from  2  to  0),  the  objective  function  (2)  with  r,  given  by  Eq.  (3) 
becomes  a  function  of  EUD.  Generally  speaking,  the  formalism  proposed  here  is  a  hybrid 
of  dose-based  and  EUD-based  functions.  Practically,  the  optimization  problem  here  is 
formulated  at  a  voxel  level  in  dose  domain  while  the  EUD-based  formalism  works  at  a 
structural  level  in  EUD  domain.  The  voxel  level  information,  such  as  non-uniform 
biology  data  in  target  and/or  sensitive  structures,  can  be  more  easily  integrated  into  the 
objective  function  proposed  here.  In  addition,  prescription  and  evaluation  of  a  treatment 
plan  in  dose  domain  other  than  in  EUD  domain  makes  the  proposed  method  clinically 
more  practical. 


Conclusion 

Inverse  planning  is  essentially  to  balance  the  competing  dosimetric  or  clinical 
requirements  of  the  individual  elements  in  a  given  system.  The  level  of  optimality  of  the 
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final  IMRT  treatment  plan  depends  not  only  on  the  performance  of  the  optimization 
algorithm  but  also,  more  fundamentally,  on  the  way  that  the  tradeoff  between  different 
competing  elements  is  defined.  While  a  voxel  in  a  patient  is  the  smallest  dosimetric 
element,  the  tradeoff  in  current  inverse  planning  algorithm  is  being  done  implicitly  on  an 
organ  or  structure  level.  In  this  work,  we  have  shown  that  this  type  of  tradeoff  scheme  is 
deficient  and  seriously  limits  the  solution  space.  Furthermore,  we  proposed  a  model  for 
the  tradeoff  of  different  voxels  through  the  modulation  of  voxel  dependent  importance 
factors.  A  novel  inverse  planning  algorithm  with  adaptively  determined  local  importance 
factors  was  described.  A  comparison  with  the  conventional  inverse  planning  technique 
indicated  that  the  new  tradeoff  strategy  substantially  improves  the  IMRT  dose 
distributions  and  allows  us  to  obtain  solutions  that  would  otherwise  be  unattainable. 
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Table  I  S  ummary  of  the  optimization  parameters  used  in  the  conventional  and  newly 
proposed  approaches  for  the  hypothetical  phantom  case 


The  dose-based  approach  The  proposed  approach 


Organs 

Importance 
factors  (ra) 

Prescription 
or  tolerance 
doses  (%) 

Importance 
factors  (ra) 

Dref(%) 

k 

Target 

8.0 

100 

15.0 

100 

-4 

Critical 

5.0 

35 

0.6 

45 

10 

Structure 
Normal  tissue 

1.0 

70 

0.5 

70 

2 

Table  II  Summary  of  the  optimization  parameters  used  in  the  conventional  and  newly 
proposed  approaches  for  the  prostate  case 


The  dose-based  approach  The  proposed  approach 


Organs 

Importance 
factors  {rG) 

Prescription 
or  tolerance 
doses  (Gy) 

Importance 
factors  (ra) 

Dref 

(Gy) 

k 

Target 

4.0 

70 

5.0 

70 

-5 

Bladder 

1.0 

50 

0.1 

65 

2 

Rectum 

1.0 

45 

0.3 

60 

8.3 

Normal  tissue 

0.5 

65 

0.4 

68 

2 
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Captions 


Fig.  1.  The  comparison  of  the  isodose  distributions  for  the  hypothetical  phantom  case:  (a) 
the  geometry  of  the  phantom;  and  (b)  dose  distribution  obtained  using  structurally 
uniform  importance  factors;  (c)  dose  distribution  obtained  using  the  newly  proposed 
algorithm  with  structurally  non-uniform  importance  factors. 

Fig.  2 .  Comparison  of  the  Dose  V olume  Histograms  (DVH)  for  IMRT  plans  obtained 
using  our  newly  proposed  algorithm  (solid  curves)  and  the  conventional  algorithm  (dash 
lines). 

Fig.  3.  Comparison  of  the  isodose  distributions  for  the  prostate  IMRT  plans:  (a)  dose 
distribution  obtained  using  structurally  uniform  importance  factors;  (b)  dose  distribution 
obtained  using  the  newly  proposed  algorithm  with  structurally  non-uniform  importance 
factors.  The  results  on  two  transverse  slices,  and  one  sagittal  slice  are  shown. 

Fig.  4.  Comparison  of  the  Dose  Volume  Histograms  (DVH)  for  the  prostate  IMRT  plans 
obtained  using  our  newly  proposed  algorithm  (solid  curves)  and  the  conventional 
algorithm  (dash  lines). 
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Abstract 


The  endorectal  coil  is  being  increasingly  used  in  magnetic  resonance  imaging  (MRI)  and 
MR  spectroscopic  imaging  (MRSI)  to  obtain  anatomic  and  metabolic  images  of  the 
prostate  with  high  signal-to-noise  ratio  (SNR).  In  practice,  however,  the  use  of  endorectal 
probe  inevitably  distorts  the  prostate  and  other  soft  tissue  organs,  making  the  analysis  and 
the  use  of  the  acquired  image  data  in  treatment  planning  difficult.  The  purpose  of  this 
work  is  to  develop  a  deformable  image  registration  algorithm  to  map  the  MRI/MRSI 
information  obtained  using  an  endorectal  probe  onto  CT  images  and  to  verify  the 
accuracy  of  the  registration  by  phantom  and  patient  studies.  A  mapping  procedure 
involved  using  a  thin  plate  spline  (TPS)  transformation  was  implemented  to  establish 
voxel-to-voxel  correspondence  between  a  reference  image  and  a  floating  image  with 
deformation.  An  elastic  phantom  with  a  number  of  implanted  fiducial  markers  was 
designed  for  the  validation  of  the  quality  of  the  registration.  Radiographic  images  of  the 
phantom  were  obtained  before  and  after  a  series  of  intentionally  introduced  distortions. 
After  mapping  the  distorted  phantom  to  the  original  one,  the  displacements  of  the 
implanted  markers  were  measured  with  respect  to  their  ideal  positions  and  the  mean  error 
was  calculated.  In  patient  studies,  CT  images  of  two  prostate  patients  were  acquired, 
followed  by  3  Tesla  (3T)  MR  images  with  a  rigid  endorectal  coil.  Registration  quality 
was  estimated  by  the  centroid  position  displacement  and  image  coincidence  index  (Cl). 
Phantom  studies  show  that  using  the  deformable  registration  method  the  mean  landmark 
displacement  error  was  0.62  ±  0.39  mm  when  the  distortion  was  of  the  order  of  23.07 
mm.  When  a  non-deformable  method  based  on  a  rigid-body  transformation  and  scaling 
was  used  for  the  same  distortion,  the  mean  displacement  of  the  fiducials  with  respect  to 
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their  actual  positions  was  found  to  be  as  large  as  12.95  ±  0.57  mm.  For  both  patient 
studies,  significantly  improved  registration  accuracy  was  also  achieved.  The  prostate 
centroid  position  displacement  was  0.58  ±0.10  mm  and  the  Cl  was  92.6  ±  5.1  %  when  a 
TPS  transformation  was  used.  Different  from  the  non-deformable  approach,  the  TPS- 
based  registration  accommodates  the  organ  distortion  and  enables  us  to  achieve 
significantly  higher  MR/MRSI  and  CT  image  registration  accuracy.  The  technique  should 
be  useful  to  map  the  MR  spectroscopic  dataset  acquired  with  ER  probe  onto  the  treatment 
planning  CT  dataset  to  guide  radiotherapy  planning. 
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1.  Introduction 


The  introduction  of  endorectal  (ER)  surface  coils  significantly  improves  the 
spatial  resolution  and  signal-to-noise  ratio  (SNR)  of  prostate  MR  and  MR  spectroscopic 
imaging  1_8.  The  new  MRI/MRSI  tool  provides  an  unprecedented  means  for  us  to 
characterize  the  location(s)  and  volume(s)  of  intraprostatic  lesion(s)  and  to  evaluate  the 
possible  capsular  penetration,  invasion  of  neurovascular  bundle,  and  seminal  vesicle 
involvement  24,  9'13.  The  information  derived  from  the  new  imaging  modality  is  also 
valuable  for  guiding  radiation  therapy  treatment  planning  to  escalate  radiation  doses 
according  to  the  regional  tumor  burden  14'18.  In  practice,  the  use  of  ER  coil  severely 
distorts  the  prostate  and  surrounding  organs.  On  the  other  hand,  current  radiation 
treatment  planning  is  performed  on  the  CT  images  without  distortion.  In  order  to  use  ER- 
based  image  data  to  guide  radiation  therapy  treatment  planning,  it  is  imperative  to 
develop  a  method  to  map  the  information  in  the  ER-based  MRI/MRSI  to  the 
corresponding  location  in  CT  images.  Zaider  et  al  (2000)  and  Mizowwki  et  al  (2002) 
have  reported  a  translation  and  scaling  based  registration  method  to  map  MRS  positive 
volumes  onto  the  CT  and  ultrasound  images.  In  their  approach,  the  coordinates  of  the 
boundary  and  the  center  of  mass  were  used  to  linearly  interpolate  the  positions  of  the 
mapped  voxels.  Although  a  clinically  acceptable  mean  difference  (2.4  mm)  between  the 
predicted  and  measured  positions  was  reported,  larger  discrepancy  was  found  for  regions 
with  more  severe  distortion  (>  4mm). 

The  purpose  of  this  paper  is  to  apply  thin  plate  spline  (TPS)  -based  deformable 
registration  to  improve  the  previously  reported  nondeformable  MRS  and  CT  mapping 
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technique  and  test  the  registration  accuracy  using  a  series  of  phantom  measurements.  The 

TPS  transformation  is  a  well-established  mathematical  method  and  its  central  idea  is  to 

find  a  continuous  transformation  to  minimize  the  difference  between  the  control  points  in 

two  images.  Since  its  first  introduction  into  medical  image  analysis  (Bookstein  1989),  the 

TPS  has  been  successfully  used  on  several  applications.  A  two  step  registration  scheme 

(rigid  body  registration  and  TPS  warping)  was  employed  to  make  comparisons  of  MR 

images  in  interventional  MRI  guided  radiofrequency  ablation  to  determine  whether  a 

tumor  is  adequately  treated  19 ’ 20 .  In  order  to  map  changes  in  the  shape  and  position  of  the 

liver  between  inhale  and  exhale  breath  held  CT  models  of  a  patient,  a  mutual  information 

(MI)  based  alignment  with  TPS  warping  was  proposed  21 .  A  TPS  transformation  based 

technique  has  also  been  found  useful  to  correct  image  distortion  in  fluoroscopic  images 
22 


2.  Methods  and  materials 

2.1.  Image  Mapping  Method 

The  prostate  volumes  in  CT  and  MRI  were  first  contoured  by  an  experienced 
oncologist.  The  rotation  operator  was  applied  to  adjust  the  relative  tilt  between  two 
volumes.  The  axial  slices  of  CT  and  MR  data  set  were  resampled  using  1  mm  interval. 
We  aligned  CT  and  MR  slices  with  reference  to  the  apex  and  base  of  the  glands.  Four  to 
eight  control  points  were  placed  in  each  pair  of  slices.  The  control  points  were  only  put 
along  the  contour  of  the  gland  and  they  are  featured  points  such  as  comers  and 
intersections  of  edges.  Lastly  the  TPS  transformation  was  applied  on  each  pair  of  slices  to 
establish  a  mapping  relationship  between  voxels  of  MRI  and  CT.  For  convenience, 
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henceforth,  the  non-distorted  CT  volume  is  referred  to  as  the  reference  and  the  distorted 
MRI  the  floating  dataset. 

The  detailed  description  of  the  TPS  transformation  can  be  found  in  Bookstein’s 
paper 23 .  Here  we  briefly  summarize  the  computational  procedure. 

(1)  Assume  Pi  =(xi,yi),  Pi  ={xi,yi),...,  P„  =  (xn,yn )  are  the  n  control  points 
in  the  reference  image.  The  distance  between  point  i  and  j  is  given  by  nj  =  \  Pi  -  Pj\. 
Define  matrices 


1  jci  yi 

p= 

1  X2  yi 

9 

1  Xn  yn 

0 

U(rn) 

U(rin) 

U(rii) 

0 

U(rm) 

U{rn  i) 

U(rm)  ••• 

0 

and 


K  P 
PT  O 


where  O  is  a  3  x3  matrix  of  zeros  and  t/  is  a  basis  function  U(r)  =  r2  log  r2 . 

(2)  Let  Q\  =(wi,vi),  Qi  ={ui,vi),  ...,  Qn  - (u„, v«)  be  n  corresponding  control 
points  in  the  floating  image.  We  construct  matrices 


V  = 


Ml  Ul  •••  Un 
VI  V2  •  •  •  Vn 


Y  =  (v\  0  0  0)r. 
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The  weighting  vector  W  =  (wi,W2,  -,wn)  and  the  coefficients  a\t  aUi  and  av  can  be 
computed  by  the  equation 

L-'Y  =  (W\  a,  a„  avf . 

(3)  Use  the  elements  of  LAY  to  define  a  function  /(m'jV1)  everywhere  in  the  plane 

f{u  ',v')  =  al  +  auu  +  avv  +  wJJ (|  Pi  -  (u,  v)  |). 

/=0 

This  function  transforms  a  voxel  in  the  floating  volume  to  a  new  coordinate  in  the 
reference  volume.  The  computation  of  the  TPS  transformation  is  rather  efficient.  In  our 
experiment,  it  took  around  5  seconds  to  compute  a  520  x  520-pixel,  8-control  point 
transformation  on  a  Personal  Computer  (PC)  with  Intel  Pentium®  II  400Mhz  CPU  (Intel 
Corporation,  Sunnyvale,  CA)  and  256MB  memory. 

For  comparison  purpose,  we  also  implemented  the  non-deformable  registration 
method  reported  by  Zaider  et  al  (2000)  and  Mizowaki  et  al  (2002).  For  a  particular  voxel 
in  the  MR  space  (coordinate  zi),  the  z  coordinate  in  the  US/CT  space  was  obtained  from: 


z\  ~  zc,  _  Z2  "  zc.2 


where  zxi  and  zT2  are  the  coordinates  of  the  superior  aspects  of  the  prostate  in  the  MR  and 
US/CT  volume,  respectively,  zBi  and  zB2  refer  to  the  z  coordinates  of  the  inferior  aspects 
of  the  prostate,  respectively,  and  zci  and  zc2  represent  the  z  coordinates  of  the  prostate 
centroid  in  the  MR  space  and  US/CT  space,  respectively.  Similarly,  the  (x,  y)  coordinates 
were  mapped  as  follows: 
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y a\  -  y*  _  y±  -  yh 


y  -  y\  y  a,  -  y2 

xl\  ~  —  *h~_  Xr2 


Here,  yA  and  yp  are  the  y  coordinates  of  the  anterior  and  posterior  aspects  of  the 
prostate,  respectively,  and  xL  and  xR  refer  to,  respectively,  the  x  coordinates  of  the  left 
and  right  aspects  of  the  prostate.  The  results  obtained  by  using  this  method  and  the  newly 
developed  TPS  method  were  quantitively  compared  in  the  phantom  and  patient  studies. 


2.2.  Phantom  Construction  and  Imaging 

Tissue  equivalent  bolus  material  was  used  to  construct  the  2D  phantoms  which 
simulate  the  axial  sections  of  the  patient  dataset.  The  bolus,  made  of  vinyl  gel,  is  elastic 
and  has  a  density  close  to  that  of  water.  Ten  to  fifteen  metal  fiducial  landmarks  were 
embedded  into  each  phantom.  The  phantoms  were  held  in  a  custom  made  plastic  holder, 
allowing  them  to  be  constricted  and  deformed  in  specifically  chosen  regions  (Fig.l).  The 
radiographic  images  were  then  acquired  in  anterior-posterior  (AP)  and  lateral  (LT) 
direction  using  a  Ximatron  Radiotherapy  Simulator  (Varian  Medical  Systems,  Palo  Alto, 
CA).  An  analysis  of  the  AP/LT  images  for  each  phantom  revealed  the  geometric 
locations  of  the  fiducial  markers  before  and  after  the  distortion. 

2.3.  Patient  Image  Acquisition 

Patient  MRI  was  acquired  on  a  3-Tesla  MR  scanner  (Signa;  GE  Medical  Systems, 
Milwaukee,  WI).  RF  excitation  was  achieved  by  using  the  whole  body  birdcage 
resonator,  and  the  MR  signal  was  received  using  a  4-element  phased-array  antenna  (G.E. 
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Medical  Systems,  Milwaukee,  WI)  combined  with  a  rigid  single  loop  receiver-only 
surface  coil  with  a  fixed  geometry  that  enables  optimal  tuning  and  matching  for  use  at 
3T.  The  coil  dimensions  are  similar  to  trans-rectal  ultrasound  transducers  used  for  routine 
sonographically-guided  prostate  imaging  and  biopsy.  The  ER-induced  distortion  of  MRSI 
is  very  close  to  that  of  MRI.  We  show  MR  image  in  this  study  because  they  have  higher 
image  quality  than  MRSI.  Patient  CT  images  were  acquired  using  a  PQ5000  CT  Scanner 
(Philips  Medical  Systems,  Cleveland,  OH). 


2.4.  Validation  of  the  Image  Registration 

For  phantom  studies,  the  control  points  were  chosen  only  in  the  periphery  for  the 
registration  of  the  floating  and  reference  images.  The  inserted  landmarks  were  used  to 
trace  the  displacement  and  verify  the  registration  accuracy.  After  mapping  the  distorted 
phantom  to  the  original  one,  the  displacements  of  the  implanted  markers  were  measured 
with  respect  to  their  ideal  positions  and  the  mean  discrepancy  was  calculated  for  each 
phantom.  The  mean  landmark  displacement  (MLD)  was  used  as  a  metric  for  evaluating 
the  quality  of  the  registration. 

For  patient  studies,  typically  6-8  control  points  were  chosen  along  the  contour  of 
the  prostate  based  on  the  pronounced  anatomical  feature.  Patient  MR  and  CT  registration 
accuracy  was  estimated  by  using  the  centroid  position  displacement  of  the  prostate  and 
the  coincidence  index  (Cl)  defined  by 


CI(R,F)  =  - 


_[F(;c)R(x)</jc 


£  max  {.F(x),  R(x)|  dx 

where  Cl  is  unity  when  two  structures  overlap  exactly  and  zero  when  they  are  completely 


disjoined  24.  The  floating  (F)  image  and  the  reference  (R)  images  were  converted  to 
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binary  for  the  calculation.  The  use  of  Cl  provided  us  with  an  effective  measure  of  the 
similarity  between  the  warped  floating  (F)  image  and  the  reference  (R)  image.  In  both 
evaluations  the  tissue  density  was  assumed  to  be  homogenous. 

3.  Results 

3.1.  Phantom  studies 

We  first  studied  the  dependence  of  registration  accuracy  on  the  number  of  control 
points.  An  elastic  phantom  with  dimension  5.5  X5.5  xl  cm3  was  used  here.  The 
phantom  was  distorted  by  the  insertion  of  an  object  in  the  holder  (Fig.  2A  left)  and  it 
restored  to  the  original  shape  when  the  object  was  removed  (Fig.  2A  right).  When  four 
control  points  were  selected  along  the  margin  (Fig.  2B  left),  we  obtained  the  warped 
image  shown  in  the  middle  column.  To  evaluate  the  TPS  algorithm,  we  computed  the 
difference  between  the  TPS  predicted  and  the  true  image  (Fig.  2B  right  panel).  The  MLD 
was  found  to  be  1.51  ±  0.49  mm.  It  is  seen  that  the  implanted  landmarks  do  not  coincide 
well  in  the  two  images.  Next  we  added  two  more  control  points  in  the  periphery  and  the 
corresponding  mapped  image  shows  reduced  registration  error  with  MLD  down  to  0.76 
±  0.54  mm  (Fig.  2C).  When  eight  control  points  were  selected,  the  MLD  was  further 
reduced  to  0.46  ±  0.34  mm  and  no  significant  landmark  displacement  was  found  in  the 
difference  image  (Fig.  2D).  In  Fig.  2E  we  summarized  the  MLDs  when  four,  six  and 
eight  control  points  were  used  in  the  warping  calculation.  The  use  of  more  control  points 
resulted  in  higher  registration  accuracy.  In  practice,  however,  increasing  the  number  of 
control  points  requires  additional  manual  interaction  and  prolongs  the  registration 
process.  In  the  following  studies,  six  to  eight  pairs  of  control  points  were  selected  for  the 
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TPS  registration.  We  also  mapped  the  distorted  phantoms  onto  the  reference  images 
using  a  rigid-body  registration  and  scaling  based  method  17,  25.  The  non-deformable 
registration  resulted  in  a  2.50  ±  0.83  mm  MLD  when  the  maximum  displacement  was 
4.2  mm.  Hence,  in  the  situation  of  a  4.2  mm  distortion,  the  TPS  method  with  eight 
control  points  yielded  an  MLD  that  was  only  1 8.4%  of  the  MLD  obtained  with  the  non- 
deformable  model  (0.46  versus  2.50  mm). 

In  the  next  level  of  validation  we  tested  the  algorithm  with  a  larger  rectangular 
phantom  (9.2 x  5.1x1  cm)  to  allow  more  flexible  distortions.  The  reference  image  is 
shown  in  Fig.  3A.  The  floating  images  under  a  few  different  levels  of  distortions  are 
shown  in  the  left  columns  of  B,  C  and  D.  Eight  control  points  were  used  here  to  register 
the  floating  and  reference  images.  The  middle  panels  of  B,  C  and  D  show  the  results  after 
the  TPS  transformation.  The  differences  between  the  TPS  predictions  and  references  are 
shown  on  the  right  panels  of  B,  C  and  D.  The  quality  of  the  TPS  mapping  was  accessed 
by  using  the  maximum  and  mean  landmark  displacement  (MLD).  As  summarized  in 
Table  1,  for  the  distortion  shown  in  Fig.  3A,  the  non-deformable  registration  gave  a  MLD 
of  4.62  ±  2.71  mm,  whereas  the  deformable  registration  reduced  the  error  down  to  0.45 
±  0.53  mm.  For  the  studies  shown  in  Fig.  3C  and  D,  the  non-deformable  registration 
yielded  MLDs  to  7.35  ±  4.20  mm  and  12.95  ±  6.57  mm,  respectively.  The  application  of 
the  deformable  warping  module  significantly  improved  the  mapping  and  led  to  MLDs  of 
0.57  ±  0.49  mm  and  0.62  ±  0.39  mm,  respectively.  The  largest  registration  error 
between  the  TPS  prediction  and  the  ideal  situation  were  found  to  be  1.09  mm,  1.05  mm 
and  0.99  mm,  respectively,  for  the  three  phantom  distortions. 
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In  order  to  examine  the  consistency  of  the  registration,  we  inversed  the  previous 
transformation  procedure  by  transforming  the  TPS-warped  images  (Fig.3  B,  C  and  D 
middle  panels)  back  to  the  distorted  floating  images.  The  calculation  results  are  shown  in 
the  left  columns  of  Fig.4  A,  B  and  C.  The  difference  images  between  them  and  the 
original  deformed  phantom  images  (the  left  columns  of  Fig.  3)  are  shown  in  the  right 
panels.  The  resultant  overlap  of  fiducial  points  was  excellent  in  all  three  cases,  suggesting 
the  TPS  is  capable  of  generating  consistent  good  mapping  independent  of  the  starting 
images.  The  MLDs  for  the  three  groups  were  0.23  ±  0.08,  0.23  ±  0.18  and  0.20  ±0.11 
mm,  respectively.  Maximum  landmark  registration  discrepancies  were  found  to  be  0.35, 
0.50  and  0.31  mm,  respectively  (Fig.  4  D). 

3.2.  Patient  studies 

We  first  studied  where  the  distortion  most  likely  happens  in  the  ER-based  MR 
images.  After  target  segmentation  and  rotation  operation,  we  compared  the  dimensions  of 
the  prostate  in  the  datasets  of  two  patients.  The  height  of  the  prostate  along  superior- 
posterior  axis  was  found  almost  the  same  (3.1%  discrepancy)  in  CT  and  ER  based-MR 
images.  The  width  along  left-right  axis  and  the  length  along  anterior-posterior  axis  differ 
a  lot  between  CT  and  MR  images.  We  measured  the  width  and  length  of  the  prostate  in 
the  middle  axial  slices  of  two  patients.  The  width  of  the  prostate  in  MRI  is  1 16.9  ±  1.5% 
of  that  of  the  prostate  in  CT.  The  length  of  the  prostate  in  MRI  is  82.1  ±  1.4%  of  that  of 
the  prostate  in  CT.  After  further  comparing  the  shape  of  the  prostates  in  CT  and  ER-MRI, 
we  conclude  that  distortion  mostly  happens  in  the  transverse  plane.  This  observation 
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suggests  that,  to  a  reasonable  approximation,  we  could  perform  the  mapping  procedure  in 
a  slice  by  slice  fashion. 

The  TPS  transformation  was  applied  to  the  coregistration  of  the  CT  and  ER-based 
MR  images.  We  show  a  representative  axial  slice  of  a  patient’s  CT  images  in  Fig.  5 A.  To 
have  a  better  view  of  the  volume  of  interest,  we  selected  a  rectangular  region 
encompassing  the  prostate  (Fig.  5 A  right).  The  MR  images  were  acquired  with  high 
resolution  and  the  posterior  portion  of  the  image  was  distorted  by  the  presence  of  the  ER 
coil  (Fig.  5B  left).  Eight  control  points  were  chosen  along  the  contour  of  the 
corresponding  MR  and  CT  images.  TPS  transformation  was  applied  to  the  distorted  MR 
image  and  the  mapped  MRI  contour  of  the  prostate  overlapped  almost  completely  with 
that  from  the  CT  scan  (Fig.  5B  right).  Difference  between  TPS-derived  MR  contour  and 
CT  contour  is  shown  in  Fig.  5C.  Most  prostate  regions  were  in  good  agreement  including 
the  seriously  contorted  left  and  right  posterior  regions  of  the  image.  Similar  results  were 
obtained  from  the  other  patient. 

We  used  centroid  position  displacement  and  coincidence  index  (Cl)  between  the 
mapped  MR  images  and  the  CT  images  to  quantify  the  registration  accuracy  for  these  two 
datasets.  Using  the  TPS  method,  the  centroid  displacement  was  0.58  ±  0.10  mm, 
significantly  less  than  that  of  the  non-deformable  registration  (1.96  ±  0.43  mm,  Fig.  6A). 
The  Cl  indices  were  found  to  be  close  to  unity  (92.6  ±  5.1  %),  indicating  that  the  TPS 
algorithm  is  able  to  model  the  nonrigid  soft  tissue  deformation  caused  by  the  endorectal 
coil  placement  (Fig.  6B).  On  the  other  hand,  a  much  lower  Cl,  50.0  ±  9.4  %,  was  found 
when  using  the  non-deformable  registration.  This  suggests  that  fusion  with  a  rigid-body 
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transformation  and  scaling  is  inadequate  to  deal  with  the  system  involving  the  images 
acquired  with  the  ER  coils. 

The  registration  error  depends  on  the  appropriate  placement  of  control  point.  We 
studied  registration  inconsistencies  between  different  trials  of  one  operator  and  between 
three  operators.  The  intra-operator  experiment  was  repeated  five  times  on  one  patient’s 
data.  The  centroid  displacement  was  found  to  be  in  the  range  of  0.31  to  0.65  mm.  The  Cl 
indices  were  found  to  be  from  91.7%  to  93.5%.  Three  operators  were  asked  to  repeat  the 
control  point  placement  five  times  and  the  mean  results  between  the  operators  were 
compared.  The  centroid  displacements  were  found  to  be  0.55  ±  0.30,  0.47  ±  0.17,  and 
0.54  ±  0.25  mm,  respectively.  The  Cl  indices,  corresponding  to  three  operators,  were 
92.7  ±  0.9%,  93.5  ±  0.7%  and  92.5  ±  1.2%,  respectively.  The  centroid  displacement 
and  Cl  index  show  no  significant  difference  between  trials  and  operators. 

4.  Discussion  and  Conclusion 

Registration  has  been  implemented  in  several  commercial  medical  image  analysis 
and  radiation  treatment  planning  systems.  For  example,  Radionics  (Radionics™, 
Burlington,  MA)  has  developed  ImageFusion  software  which  provides  the  ability  to  fuse 
multiple  image  sets  based  on  the  mutual  information.  AcQSim  Oncodiagnostic 
Simulation/Localization  System  (Philips  Medical  Systems,  Cleveland,  OH)  provides  two 
registration  methods:  point  matching  (a  minimum  of  three  common  points  need  to  be 
selected  on  both  sets  of  images  registered)  and  interactive  image-based  registration  (a 
color  wash  of  one  image  set  is  displayed  over  a  grayscale  image  of  the  other).  At  this 
point,  they  all  use  a  rigid-body  transformation  and  scaling,  which  maintain  the 
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straightness  of  lines,  and  hence  cannot  accommodate  contour/shape  distortion.  In  reality, 
the  shape  of  the  prostate  gland  can  be  easily  changed  by  many  factors  such  as  patient 
position  change,  invasive  brachytherapy  procedures  or  endorectal  coil  placement  during 
high  resolution  MR/MRS  images  acquisition  26.  To  help  physicians  to  segment  the 
prostate  gland  and  possible  intraprostatic  lesions  by  incorporating  MRI/MRSI  metabolic 
data  on  a  CT-based  treatment  planning  system,  there  is  an  indisputable  need  for 
developing  a  computationally  efficient  deformable  registration  technique  to  achieve 
voxel  to  voxel  mapping.  In  this  work,  we  used  a  TPS  method  to  register  the  endorectal 
coil-based  MR  data  with  CT  images.  The  data  presented  in  the  last  section  suggests  that 
the  TPS  technique  is  well  suited  for  this  type  of  application. 

The  warping  process  was  carried  out  in  a  slice-by-slice  way  and  is  worth  of 
further  investigation.  This  may  result  in  the  registration  error  in  the  longitudinal  direction. 
Based  on  our  observations  for  the  two  patients  involved  in  this  study,  it  seems  that  the 
distortion  occurs  mainly  along  the  right-left  and  the  anterior-posterior  directions.  The 
height  of  the  prostate  along  the  superior-inferior  axis  remains  almost  constant  in  the  MR 
and  CT  datasets.  This  is  consistent  to  the  finding  by  another  group  in  1 .5T  MR  imaging 
of  the  prostate  26 .  In  actuality,  it  is  possible  to  extend  the  current  quasi-2D  model  to  a 
fully  3D  one.  The  current  study  sheds  useful  insight  into  this  type  of  extension  and 
provides  a  nature  starting  point  for  the  implementation  of  a  complete  3D  TPS  mapping. 

A  few  more  sophisticated  deformable  registration  methods  have  been  investigated 
by  several  groups.  A  viscous-  fluid  transformation  and  fluid-  landmark  registration 
technique  have  been  proposed  to  model  the  nonrigid  deformation  of  organs  in 
intracavitary  brachytherapy  24,  27 .  A  finite-element  method  has  been  used  to  model  the 
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tissue  mechanical  property  and  to  register  brain  and  prostate  images  28,  29 .  A 
biomechanical  model  of  an  elastic  body  has  been  used  to  quantify  patient  organ  motion  in 
the  process  of  radiation  therapy  so  that  the  dose  delivered  on  the  volume  of  a  deforming 
organ  can  be  accumulated  30.  These  methods  are  usually  computationally  intensive. 
Moreover,  the  model  parameters  normally  need  to  be  determined  empirically  because  of 
the  lack  of  tissue  biomechanical  data  in  the  literature,  which  compromises  the  advantages 
of  these  physics-based  models.  The  TPS  method  is  a  simple  deformable  registration 
technique  based  on  the  mathematical  interpolation.  Our  phantom  and  patient  studies  have 
shown  that  the  approach  is  computationally  efficient  and  can  yield  clinically  acceptable 
registration  accuracy. 

It’s  noted  that  the  TPS  based  registration  needs  manual  placement  of  control 
points,  which  requires  the  input  from  an  experienced  clinician.  This  is  similar  to  the 
previously  reported  rigid  body-  based  registration  method  17,  2S.  The  intra-  and  inter¬ 
operator  experiments  have  shown  that  the  registration  accuracy  doesn’t  not  depend  on  the 
different  operators  or  different  trials  significantly. 

In  conclusion,  we  have  implemented  a  TPS  transformation  algorithm  to  map 
voxels  in  endorectal  coil-based  prostate  MR/MRS  images  with  those  in  CT  images.  The 
deformable  mapping  technique  significantly  improved  the  previously  reported  non- 
deformable  method  and  should  be  adequate  for  routine  clinical  application.  The  accuracy 
of  the  approach  has  been  tested  by  using  phantom  and  patient  studies.  The  registration 
scheme  should  be  useful  to  map  the  functional  MRSI  data  onto  CT  to  guide  the  design  of 
conformal  radiation  treatment  plans. 
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Table  1 


Deformable  and  non-deformable  registration  error  of  the  distorted  phantoms  shown  in 
Fig.  3  B,  C  and  D.  The  maximum  refers  to  the  maximum  landmark  displacement  and  the 
mean  refers  to  the  mean  landmark  displacement  in  each  case. 


Distortion 

(mm) 

Non-deformable  registration 

Deformable  registration 

maximum 

mean 

maximum 

mean 

Fig.  3B 

9.83 

4.62  ±  2.71 

1.09 

0.45  ±  0.53 

Fig.  3C 

14.74 

7.35  ±  4.20 

1.05 

0.57  ±  0.49 

Fig.  3D 

23.07 

12.95  ±  6.57 

0.99 

0.62  ±  0.39 
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Figure  captions 


Figure  1.  A  photo  of  the  deformable  phantom  with  implanted  landmarks  in  a  holder. 

Figure  2.  A  registration  study  using  a  square  phantom  deformed  by  external  force.  (A) 
The  phantom  under  the  influence  of  a  force  (left)  and  its  original  shape  (right).  The 
distorted  phantom  is  shown  in  a  smaller  scale  than  the  original  phantom  in  order  to 
include  part  of  the  holder  and  external  object.  (B)  The  position  of  four  control  points  on 
the  distorted  phantom  as  indicated  by  pink  plus  signs  (left).  The  middle  and  right  show 
the  computed  deformed  image  and  difference  image,  respectively.  (C)  and  (D)  are  similar 
to  B  except  that  six  and  eight  control  points  are  used  respectively.  (E)  The  landmark 
displacement  of  the  three  groups. 

Figure  3.  A  registration  study  by  systematically  bending  a  rectangular  phantom.  (A)  The 
original  phantom.  (B)  The  distorted  phantom  (left),  the  TPS-warped  phantom  (middle) 
and  the  difference  image  (right).  (C)  and  (D)  are  similar  to  (B)  except  that  with  increased 
distortions. 

Figure  4.  Registration  consistency  test.  Left  panels  of  A,  B  and  C  represent  the  computer- 
warped  images  with  the  middle  panel  images  of  Fig.  3  B,  C  and  D  as  input.  Right  panels 
of  A,  B  and  C  represent  the  corresponding  difference  images  between  the  mapped  and  the 
original  images.  D,  Landmark  displacement  between  the  model  prediction  and  the  actual 
position  for  the  three  groups. 


22 


Figure  5.  Deformable  registration  of  the  prostate  gland  in  a  patient.  (A)  A  transverse  CT 
study  (left)  and  a  region  of  interest  encompassing  the  prostate  (right).  (B)  The  MRI  study 
(left)  and  computer-deformed  image  (right).  (C)  Difference  between  the  CT  and  the 
mapped  MR  image.  The  control  points  are  denoted  with  plus  signs. 

Figure  6.  Centroid  position  displacement  (A)  and  coincidence  index  of  deformable  and 
non-deformable  registration  in  patient  studies  (B). 
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Purpose:  IMRT  is  being  increasingly  used  for  treatment  of  prostate  cancer.  In  practice,  however, 
the  beam  orientations  used  for  the  treatments  are  still  selected  empirically  without  any  guideline. 
The  purpose  of  this  work  was  to  investigate  inter-patient  variation  of  the  optimal  beam 
configuration  and  to  facilitate  IMRT  prostate  treatment  planning  by  deriving  a  set  of  beam 
orientation  class-solutions  for  a  range  of  numbers  of  incident  beams. 

Method:  Fifteen  prostate  cases  were  used  to  generate  the  beam  orientation  class-solutions.  For 
each  patient  and  a  given  number  of  incident  beams,  a  multiobjective  optimization  engine,  which 
optimizes  the  beam  directions  by  using  a  genetic  algorithm  and  the  beam  profiles  by  a  gradient- 
based  algorithm,  was  employed  to  provide  the  patient  specific  beam  directions.  For  each  pre¬ 
selected  number  of  beams,  the  angular  distributions  of  the  incident  beams  of  the  fifteen  patients 
were  analyzed  and  a  population-based  average  of  the  angular  distributions  was  conducted.  The 
most  frequent  beam  directions  are  identified  as  the  class-solution  for  the  corresponding  number  of 
incident  beams.  The  level  of  validity  of  the  class-solutions  was  tested  using  an  additional  clinical 
prostate  case  by  comparing  with  the  individually  optimized  beam  configurations. 

Results:  In  the  case  of  five  incident  beams,  beam  configuration  with  gantry  angles  of  55°,  120°, 
195°,  270°  and  345°  seems  to  yield  the  best  possible  IMRT  dose  distribution.  For  the  fifteen  cases 
considered,  the  gantry  angle  of  any  of  the  five  beams  was  all  distributed  within  a  range  of  18°, 
indicating  that  a  population-based  class-solution  may  provide  a  good  approximation  of  the  optimal 
beam  orientations.  In  other  cases  with  six  to  nine  incident  beams,  similar  trend  was  observed.  The 
class-solution  gantry  angles  for  prostate  IMRT  were  found  to  be:  3  beams-(110°,  230°,  350°),  5 
beams-(50°,  120°,  190°,  260°,  330°),  6  beams-(60°,  120°,  180°,  240°,  300°,  360°),  7  beams-(20°, 
70°,  120°,  170°,  220°,  270°,  320°),  8  beams-(20°,  65°,  110°,  155°,  200°,  245°,  290°,335°),  and  9 
beams-(30°,  70°,  110°,  150°,  190°,  230°,  270°,  310°,  350°).  The  level  of  accuracy  of  the  class- 
solution  approach  was  examined  by  analyzing  the  difference  between  the  plans  obtained  with 
class-solutions  and  patient  specific  optimizations. 

Conclusion:  A  set  of  beam  orientation  class-solutions  for  prostate  IMRT  treatment  planning  are 
provided  based  on  a  beam  orientation  optimization  study.  This  simplifies  the  trial-and-error 
selection  of  the  beam  orientations  and  facilitates  the  IMRT  prostate  treatment  planning  process. 


Key  Words:  IMRT,  class  solution,  prostate  cancer,  inverse  treatment  planning,  dose 
optimization,  beam  orientation. 
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Introduction 


Beam  configuration  may  have  significant  influence  on  the  quality  of  an  IMRT  treatment  even 
when  a  large  number  of  incident  beams  (e.g.,  9  beams)  are  used  [1].  Clinically,  however,  the 
gantry  angles  are  still  selected  empirically  and  multiple  adjustments  may  often  be  required.  This 
has  become  one  of  the  bottleneck  problems  that  hinder  the  optimal  and  efficient  use  of  IMRT.  A 
class-solution  method  for  beam  placement  has  been  thought  [2]  [3].  The  idea  behind  this 
approach  is  to  construct  a  representative  beam  configuration  based  on  previous  experience,  or 
more  adequately,  on  accurate  beam  orientation  optimization,  for  a  given  disease  site  and  then  use 
this  “class-solution”  for  subsequent  patient  treatment  planning.  Unfortunately,  up  to  this  point,  the 
degree  of  validity  of  this  approach  has  never  been  tested,  even  for  geometrically  relatively  regular 
disease  site  such  as  prostate,  where  it  is  believed  that  the  class-solution  should  work.  This  paper 
systematically  investigates  the  issue  and  derives  a  set  of  beam  orientation  class-solutions  for 
IMRT  prostate  irradiation. 

In  the  past  few  years,  several  groups  have  investigated  the  problem  of  beam  orientations  in 
3D  conformal  radiation  therapy  [4-1 1]  and  in  IMRT  [12-18].  The  main  focus  of  these  studies  was 
on  the  modeling,  formulation  and  optimization  of  beam  configurations.  To  obtain  a  set  of  good 
gantry  angles  for  an  IMRT  treatment,  in  principle,  all  one  needs  to  do  is  to  add  the  gantry  angle 
variables  into  an  objective  function  and  then  to  optimize  the  objective  function  with  respect  to  the 
gantiy  angles  and  the  beamlet  weights.  A  main  difficulty  of  this  type  of  brute-force  optimization 
lies  in  its  excessive  computational  time  caused  by  the  greatly  enlarged  search  space.  The 
development  of  beam  orientation  class-solution  for  some  simple  disease  sites  like  prostate  has  the 
potential  to  significantly  reduce  the  effort  of  beam  placement  in  IMRT  planning,  which  is  thus  of 
practical  importance. 

To  derive  a  population-based  beam  orientation  class-solution,  the  optimal  solutions  need  to 
be  obtained  for  at  least  a  representative  group  of  patients.  For  this  purpose,  in  the  next  section  we 
will  present  an  effective  beam  orientation  algorithm  for  the  optimization  of  IMRT  beam 
orientations.  It  is  intuitively  conceivable  that  the  optimal  beam  configuration  varies  from  patient 
to  patient  because  of  their  geometric  differences.  Generally  speaking,  a  class-solution  exists  if  the 
inter-patient  variation  in  the  beam  orientations  distribution  is  small.  As  will  be  seen  in  the  Results 
section,  this  is  approximately  the  case  for  the  representative  group  of  fifteen  prostate  cases 
selected  randomly  from  the  patient  pool  treated  at  Stanford.  We  have  investigated  the  behavior  of 
the  system  for  a  variety  of  numbers  of  incident  beams  and  provided  the  class-solutions  for  the 
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cases  of  five  to  nine  incident  beams  to  meet  different  clinical  needs.  Our  study  has  indicated  that 
the  beam  orientation  class-solutions  provided  here  represent  good  compromise  between 
practicality  and  optimality  and  may  facilitate  the  IMRT  treatment  planning  by  eliminating  or 
simplifying  the  trial-and-error  process  of  the  beam  orientation  selection. 

Materials  and  methods 
Multiobjective  optimization 

A  multi-objective  (MO)  optimization  with  a  mix  of  genetic  and  iterative  algorithms  was  used  to 
optimize  the  system.  The  genetic  algorithm  optimizes  the  beam  orientations  and  the  gradient 
based  method  the  intensity  profiles  of  the  beams.  The  MO  is  useful  in  radiotherapy  dose 
optimization  because  of  the  existence  of  multiple  conflicting  objectives  in  the  system.  In  the 
previous  studies,  the  problem  has  been  transformed  into  a  single  objective  problem  by  introducing 
a  set  of  importance  factors  to  parameterize  our  tradeoff  strategy  [19].  The  MO  optimization  tries 
to  obtain  all  acceptable  solutions  and  the  more  information  on  the  tradeoff  between  the  different 
objectives.  Here,  by  an  acceptable  solution  we  mean  a  plan  with  a  good  compromise  of  all  the 
objectives  involved  in  the  problem.  Mathematically,  the  MO  optimization  or  vector  optimization 
[20]  is  to  determine  a  set  of  decision  variables  that  optimizes  a  vector  function  whose  elements 
represent  M  objective  functions  without  violating  the  system  constraints. 

We  denote  the  decision  variables  in  an  optimization  problem  by  x={Xj,  j=l,2,...^V}.  To 
access  the  goodness  of  a  solution  it  is  needed  to  establish  some  criteria  of  evaluation.  These 
criteria  are  expressed  as  computable  functions  fi(x),...,  fM(x)  of  the  decision  variables,  which  are 
called  objective  functions.  These  form  a  vector  function  f.  The  vector  function  f(x)  maps  the  set  X 
to  the  set  F  that  represents  all  possible  values  of  the  objective  functions.  In  general,  some  of  these 
objectives  will  be  in  conflict  with  others  and  do  not  have  a  situation  in  which  all  the  f(x)  reach 
their  optimums  at  a  common  point  x.  The  multiobjective  optimization  problem  is  to  find  the 
vector  x=(xi,x2,...,xN),  i.e.,  the  solution  which  optimize  the  vector  function  f. 

To  proceed,  it  is  necessary  to  establish  some  criteria  to  define  the  optimality.  These 
solutions  in  multiobjective  optimization  constitute  the  so-called  Pareto  optimum  [21].  A  solution 
Xi  dominates  a  solution  x2  if  the  two  following  conditions  are  true: 

xi  is  no  worse  than  x2  in  all  objectives,  i.e.  /y(x,)<  /)(x2)  Vj  =  1 . M 

xi  is  strictly  better  than  x2  in  at  least  one  objective,  i.e.  /,(x,)<  /,(x2)  for  at  least  one 
j  e  {1,2 . M}. 
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We  assume,  without  loss  of  generality,  that  the  problem  that  we  are  dealing  with  here  is  a 
minimization  problem.  Among  a  set  of  solutions  P,  the  non-dominated  set  of  solutions  P’  are  those 
that  are  not  dominated  by  any  other  member  of  the  set  P.  When  the  set  P  is  the  entire  feasible 
search  space  then  the  set  P’  is  called  the  global  Pareto  optimal  set.  If  for  every  member  x  of  a  set 
P  there  exists  no  solution  in  the  neighborhood  of  x  then  the  solutions  of  P  form  a  local  Pareto 
optimal  set.  The  image  of  the  Pareto  optimal  set  is  called  the  Pareto  front. 

Optimization  method 

The  variables  to  be  optimized  include  the  beam  directions  and  the  bixel  intensities.  We  quantify 
the  variables  in  a  two-component  chromosome.  The  first  part  W  contains  the  bixel  weights  of  all 
beams  with  a  double  precision  floating-point  representation.  The  second  part  B  is  a  boolean  string, 
which  represents  the  possible  beam  directions.  The  boolean  value  records  if  the  beam  is  active  or 
not.  In  this  study,  the  beams  of  coplanar  setup  are  divided  in  360  beams  with  1°  separation 
between  adjacent  beams.  Each  beam  is  divided  in  a  grid  of  40  x  40  bixels  and  each  bixel  has  a 
dimension  of  lxl  cm2.  Only  bixels  that  pass  through  the  PTV  are  activated  and  considered  in  the 
optimization  process  and  the  rest  bixels  are  assigned  with  zero  weights.  The  number  of  bixels  to 
be  optimized  depends  on  the  patient’s  geometry  and  number  of  beams  ranging  from  4050  to  9584 
for  the  prostate  cases  analyzed  in  this  paper. 

The  system  was  optimized  using  two  algorithms:  a  genetic  algorithm  as  the  main 
optimization  framework  and  a  gradient-based  algorithm  applied  on  a  fraction  of  the  newly 
generated  solutions.  Without  the  gradient-based  algorithm,  the  genetic  algorithm  will  converge 
very  slowly  because  the  cost  function  is  too  sensitive  to  the  changes  in  bixel  weights  generated  by 
the  genetic  algorithm  [22].  The  gradient-based  algorithm,  used  for  a  few  iterations  on  the  new 
solutions,  lower  the  sensitivity  by  adjusting  the  bixel  weights.  While  it  is  possible  to  use  gradient- 
based  algorithms  alone  [23,  24],  a  combination  of  a  genetic  algorithm  with  a  gradient-based 
algorithm  improves  the  convergence  behavior  and  speeds  up  the  optimization.  The  combinational 
optimization  algorithm  starts  each  time  from  an  almost  optimized  set  of  bixels  taken  from  the 
parent  population,  rather  than  from  a  new  randomly  selected  set. 

Objective  functions 

The  objective  functions  used  in  this  study  for  the  PTV  and  for  an  OAR  are  defined  as 

1  Npw  / 
iv  PTV  J=\ 
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p.,  x  f)P7V  dpw 

where  ^  W  is  the  Heaviside  step  function,  Kf  is  the  prescription  dose  to  the  PTV,  °J  and 

d0AR  n 

J  are  the  calculated  doses  at  the  jth  sampling  point  for  the  PTV  and  OAR,  respectively,  NpTV 

and  ^0AR  are  the  corresponding  number  of  sampling  points.  We  include  as  an  additional  objective 

the  number  of  beams  . 

All  calculations  in  this  paper  are  performed  in  3D  geometry  with  15  MV  photons.  The 
incident  photon  fluence  was  divided  into  pencils  and  assigned  with  different  transmission  factors. 
The  width  of  the  pencil  projected  at  the  isocenter  plane  was  10  mm.  The  dose  value  at  each 
sampling  point  was  obtained  from  the  interpolation  of  tissue  maximum  ratio  (TMR)  values  and 
off-axis  ratio  (OAR)  values.  The  sampling  points  where  generated  quasi-random  in  each  structure 
[25],  their  number  corresponding  to  the  structures  volume,  at  a  density  of  30  points/cm3.  A 
minimum  number  of  1000  and  a  maximum  number  of  30000  where  the  limits  imposed 
additionally  on  the  number  of  points,  in  order  to  avoid  over-sampling  of  the  normal  tissue  or 
under-sampling  of  small  structures. 

Dosimetric  Constraints 

In  accordance  to  the  clinical  practice,  we  impose  two  sets  of  constraints  during  the  optimization 
process.  The  first  restricts  the  maximum  dose  to  an  OAR  to  a  value  below  DUmll .  The  second  one 
is  fractional-volume  limit,  since  for  some  OARs  late  effects  would  be  resulted  if  a  volume  larger 
than  a  specific  limit  VCrWc  is  irradiated  to  a  dose  exceeding  a  critical  value  DCrWc  The  dose  limits 
for  the  rectum,  bladder  and  femurs  are  assigned  according  to  the  literature  [26,  27].  For  the 
OARs,  the  selected  values  of  DUmll  and  VCrilic,DCrilic  are  summarized  in  Table  1.  The  constraints 

imposed  on  the  PTV  include  (a)  the  maximum  dose  should  be  less  than  78  cGy;  and  (b)  only  5% 
of  the  PTV  volume  may  receive  a  dose  of  less  than  72  cGy. 

The  constraints  are  implemented  by  assigning  a  better  rank  to  the  solutions  satisfying  the 
constraints,  forcing  selection  of  clinically  sensible  solutions.  This  was  realized  by  using  the 
constrained  non-domination  relation,  which  is  used  for  the  selection  of  the  individuals  of  the 
population.  The  population  in  the  NSGA-IIC  algorithm  is  sorted  according  to  the  dominance. 
Non-dominated  members  are  assigned  rank  1.  The  remaining  population  is  then  again  checked 
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and  the  non-dominated  are  assigned  rank  2.  This  process  continues  until  all  members  are 
classified. 


Optimization  calculation 

We  use  for  inverse  planning  the  MO  multi-objective  evolutionary  algorithm  non-dominated 
sorting  genetic  algorithm  with  constrained  elitism  NSGA-IIC  [28]  that  is  one  of  the  most  effective 
evolutionaiy  optimization  algorithms.  The  gradient-based  algorithm  selected  was  the  low  memory 
BFGS  algorithm  (L-BFGS)  to  reduce  the  memory  requirements  in  the  iterative  calculation  of  the 
cost  functions  gradient  [29] 

The  crossover  and  mutation  probability  Pc  and  Pm  were  chosen  as  0.9  and  0.01  for  this 
study.  The  algorithm  uses  the  arithmetic  crossover  and  the  non-uniform  mutation.  Because  this  is 
a  N-dimensional  problem  where  N  is  large,  we  use  a  population  size  of  100  to  maintain  the 
genetic  algorithm’s  diversity.  These  optimal  parameters  have  been  estimated  by  comparing  with 
results  obtained  using  a  deterministic  gradient-based  algorithm,  which  indicated  that  the  Pareto 
front  does  not  change  significantly  after  100  generations.  The  deterministic  gradient  based 
algorithm  LBFGS  is  applied  on  all  new  solutions  at  each  generation,  and  is  used  for  10  iterations. 
Without  support  from  LBFGS,  the  algorithm  would  require  1000  or  more  generations  to  produce 
similar  results. 

Selection  of  solutions  from  the  Pareto  set 

A  set  of  solutions  is  generated  at  the  end  of  a  multiobjective  optimization,  each  corresponding  to 
an  IMRT  plan  that  is  relevant  to  the  clinical  treatment  [14,21].  Specifically,  the  treatment  planner 
is  provided  with  a  table  of  values  for  all  the  members  of  the  Pareto  set  of  the  objectives,  DVHs  for 
all  OARs  and  PTV,  and  dose  statistics.  For  a  given  patient  and  pre-selected  number  of  beams,  an 
appropriate  solution  is  generally  chosen  from  the  Pareto  set  based  on  clinical  considerations  [14, 
21].  The  optimization  objectives  were  to  conform  the  PTV,  while  minimizes  the  dose  to  bladder 
and  femoral  heads. 

Beam  orientation  class-solutions 

Fifteen  prostate  cases  were  selected  randomly  from  the  patient  pool  at  Stanford  University 
Hospital.  The  patients  were  scanned  in  supine  position,  using  a  3  or  4  mm  CT  slice  thickness.  The 
contours  were  delineated  on  a  Philips  virtual  simulation  system  (Philips  Medical  System, 
Cleveland,  OH)  and  imported  in  our  in-house  treatment  planning  system  using  the  DICOM 
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transfer  protocol.  The  beam  orientation  optimization  was  then  performed  for  each  of  the  patients 
for  a  range  of  number  of  incident  beams  (N=3,  5,  6,  7,  8  and  9)  using  the  algorithm  described 
above.  Only  coplanar  beams  were  considered  in  this  work.  The  angular  distribution  of  the  incident 
beams  for  a  given  number  of  incident  beams  for  the  fifteen  patients  was  analyzed.  For  each 
number  of  beams,  we  identified  the  directions  for  the  class-solution  by  averaging  the  gantry  angles 
of  the  fifteen  patients. 

Test  of  the  degree  of  validity  of  class-solution 

In  order  to  test  the  validity  of  the  beam  orientation  class-solution  approach  we  randomly  selected 
additional  prostate  patient  from  the  patient  pool  at  Stanford  University  Hospital.  Two  IMRT  plans 
were  generated  for  the  patient:  one  with  class-solution  selected  beams  and  one  with  beams  chosen 
by  the  patient  specific  beam  orientation  optimization.  The  difference  between  the  two  plans  is 
analyzed  according  to  the  commonly  used  plan  evaluation  indices,  such  as  DVH, 
maximum/minimum/average  target/sensitive  structure  doses. 

Results 

Pareto  front  plans  for  a  given  patient 

For  a  given  patient,  a  set  of  plans  is  provided  by  the  multiobjective  optimization  described  in  the 
last  section,  that  is,  the  user  is  provided  with  a  spectrum  of  plan  index  values,  each  characterizes 
the  quality  of  a  plan  in  the  Pareto  front.  The  DVH  of  a  given  structure  also  constitutes  a  spectrum 
of  curves  instead  of  a  single  curve.  In  Figure  1  we  show  the  DVH  spectra  of  prostate,  rectum,  and 
bladder  for  a  representative  patient  for  all  the  non-dominated  five-field  IMRT  plans  in  the  Pareto 
front. 

Each  non-dominated  plan  in  the  Pareto  front  has  different  beam  parameters  (bixel  maps 
and/or  gantiy  angles).  The  angular  distributions  of  the  Pareto  front  plans  for  the  above  patient  are 
shown  in  Fig.  2  for  different  number  of  incident  beams.  If  a  beam  direction  is  more  feasible  it 
shows  up  in  the  final  plans  with  high  frequency.  For  3  beams  the  directions  appeared  with  high 
frequencies  in  the  Pareto  front  plans  are  at  120°,  240°  and  340°.  It  is  interesting  to  note  that  the 
beam  configuration  is  similar  to  that  found  by  inspecting  the  most-selected  plan.  Small  peaks 
appear  at  60°  and  290°  because  they  are  opposite  to  the  most  favorable  beam  directions,  i.e.,  240° 
and  120°. 

As  will  be  seen  in  the  next  section,  for  3-7  beams,  the  angles  of  the  maximums  of  the 
angular  distributions  are  in  accordance  with  directions  found  using  the  most-selected  plans.  As  the 
number  of  beams  increased,  the  importance  of  optimizing  beam  orientation  decreases  and  a  large 
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set  of  beam  configurations  with  more  “diffusive”  angular  distributions  become  capable  of 
producing  valid  plans  with  adequate  PTV  coverage  and  OAR  sparing.  In  this  situation,  the  feature 
of  distinct  peaks  in  the  angular  distribution  diminishes. 

Class-solutions 

For  a  given  patient,  the  beam  configuration  occurring  most  frequently  in  the  Pareto  front  is 
identified  as  the  optimal  one.  In  figure  3  we  show  the  optimized  five-beam  configurations  for  each 
of  the  fifteen  patients.  Two  symmetrical  configurations  are  observed,  with  one  beam  passing 
through  either  left  or  right  femur.  To  better  visualize  this,  in  figure  4b  we  show  the  gantry  angles 
for  the  fifteen  patients  (Figure  4)  on  a  single  plot.  We  select  the  configuration  with  one  beam 
passing  through  the  right  femur.  The  gantry  angles  of  the  ideal  class-solution  for  IMRT  treatment 
with  five  incident  beams  seems  to  be  55°,  120°,  195°,  270°  and  345°  and  are  indicated  in  red  in 
figure  4b.  These  angles  were  selected  based  on  the  average  direction  for  each  cluster  of  beam 
orientations  in  the  most-selected  plans.  For  practice  purpose,  it  is  desirable  to  “round”  the  beam 
angles  to  a  set  of  more  convenient  implementable  beams.  With  the  consideration  of  angular 
distribution  of  each  beam  direction  shown  in  figure  4b,  we  identify  5-beam  class-solution  to  be 
50°,  120°,  190°,  260°,  and  330°.  Because  of  the  adjustment  is  done  with  the  range  of  the  angular 
distribution  of  the  beams,  the  compromise  is  limited.  We  note  that  the  equispaced  setup  is  in 
accordance  to  the  findings  of  [17,  30],  and  the  starting  gantry  configuration  of  60°  is  consistent 
with  that  of  Ref.  [31]. 

The  solutions  derived  above  seem  to  be  physically  sensible.  Setting  the  angles  of  beams  1  2 
and  4  to  50°,  120°  and  260°  balances  the  dose  to  the  femoral  heads.  The  locations  of  the  3rd  (190°) 
and  5th  (330°)  beams  are  chosen  to  balance  the  dose  requirements  of  the  PTV  and  rectum. 
Depending  on  the  geometry  of  the  case,  the  optimizer  makes  individual  adjustments  of  the  beam 
angles,  but  the  selected  directions  are  within  a  ±18°  window.  The  standard  deviations  of  the  five 
fields  are  10°,  30°,  30°,  35°,  and  25°,  respectively. 

The  same  calculation  described  above  is  repeated  for  3,  6,  7,  8  and  9  incident  beams  and  the 
distributions  of  the  optimal  configurations  for  the  fifteen  cases  are  shown  in  figure  4.  The  class- 
solution  for  each  of  the  situations  is  presented  in  figure  4  in  red.  For  convenience,  the  beam  angles 
of  the  class-solution  are  summarized  here:  3  beams-(110°,  230°,  350°),  5  beams-(50°,  120°,  190°, 
260°,  330°),  6  beams-(60°,  120°,  180°,  240°,  300°,  360°),  7  beams-(20°,  70°,  120°,  170°,  220°,  270°, 
320°),  8  beams-(20°,  65°,  110°,  155°,  200°,  245°,  290°,335°),  and  9  beams-(30°,  70°,  110°,  150°, 
190°,  230°,  270°,  310°,  350°).  After  adjustments  with  consideration  of  the  3-7  beams,  the  locations 
of  the  maximums  of  the  angular  distributions  for  the  fifteen  patients  are  consistent  with  the 
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directions  found  using  the  most-selected  plans.  The  beam  setups  tend  to  be  equispaced  when  3-7 
beams  are  used,  with  the  first  beam  placed  around  20-50°.  We  noticed  that  as  the  number  of 
incident  beams  increases,  the  beams  tend  to  be  less  confined  to  a  few  fixed  directions  (see  figure 
4).  Instead,  the  optimal  configuration  for  a  case  may  differ  from  another  significantly.  In  addition, 
the  quality  of  the  final  plan  becomes  less  sensitive  to  the  beam  angle  variations.  In  other  words,  as 
the  number  of  beams  increases,  the  importance  of  optimizing  beam  orientation  decreases.  As  a 
result,  it  is  less  necessary  to  derive  a  class-solution  for  IMRT  prostate  treatment  with  8  or  9 
beams.  Practically,  however,  it  is  still  desirable  to  have  a  set  of  fixed  gantry  angles  to  “streamline” 
the  clinical  IMRT  treatment  planning  process.  For  this  reason,  we  have  specified  a  set  of  beam 
configurations  for  these  situations  in  Figure  4. 

It  is  interesting  to  observe  that,  in  the  case  of  6  and  8  beams,  one  pair  of  beams  are 
required  to  incident  from  opposed  (or  almost  opposed)  directions  in  the  class-solution.  We  expect 
that  when  lower  energy  photons  are  used,  more  opposed  beam  pairs  will  show  up  in  the  final 
solution. 

Evaluation  of  class-solution  approach 

Validity  of  the  class-solution  approach  was  invertigated  by  comparing  plans  obtained  based  on 
patient  specific  beam  orientation  optimization  and  beam  orientation  class-solution.  An  additional 
patient  was  selected  from  the  patient  pool  of  Stanford  University  Hospital.  Two  sets  of  plans  were 
generated.  In  the  first  set,  the  optimizer  uses  fixed  beams  during  the  optimizatiton  proces,  with  the 
directions  taken  from  the  class-solution,  while  in  the  second  set  beams  orientation  are  optimized. 
Both  5  and  9  incident  beams  were  studied. 

In  figure  5a  we  show  the  Pareto  fronts  of  the  two  plans  when  5  beams  are  used.  As  can  be 
seen  from  the  figure,  the  Pareto  fronts  obtained  with  the  two  methods  are  similar.  For  comparison, 
we  selected  one  optimal  plan  from  each  set  and  made  a  side-by-side  compariosn  (Figure  5b).  The 
two  plans  were  found  to  be  quite  similar  and  both  plans  satisfied  the  imposed  constraints.  The  plan 
with  class-solution  beam  configuration  irradiates  slightly  more  volume  of  rectum  at  lower  doses 
but  is  slighly  better  at  high  doses.  Similar  results  were  found  in  the  case  of  9  incident  beams,  as 
can  be  seen  from  figure  6. 

Conclusion 

In  general,  the  beam  orientation  class-solution  is  an  oversimplified  approach  since, 
clinically,  the  anatomy,  target  shape  and  its  relative  location  with  respect  to  the  adjacent  sensitive 
structures,  and  the  patient's  radition  treatment  history  may  vary  from  patient  to  patient.  An 
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individualized  beam  configuration  is  frequently,  if  not  always,  needed  to  achieve  the  best  possible 
treatment.  In  actuality,  however,  there  exist  a  few  disease  sites  (such  as  prostate)  where  the 
shapes,  locations  and  prescribed  doses  of  various  structures  involved  in  the  treatment  vary  slightly 
from  case  to  case.  In  this  situation,  beam  angle  class-solution  approach  may  be  an  acceptable 
compromise.  In  this  work,  we  have  studied  the  inter-patient  variation  of  the  optimal  beam 
configurations  for  prostate  IMRT  and  derived  beam  orientation  class-solutions  for  the  cases  of  5-9 
incident  fields  based  on  analysis  of  patient  specific  beam  orientation  optimization  for  fifteen 
prostate  cases.  The  level  of  accuracy  of  the  class-solution  approach  was  examined  by  analyzing 
the  difference  between  the  plans  obtained  with  class-solutions  and  patient  specific  optimizations. 
We  found  that,  when  the  number  of  incident  beams  is  less  or  equal  to  7,  the  angle  variations  from 
patient  to  patient  are  within  18,  suggesting  the  existence  of  beam  angle  class-solution  for  prostate 
IMRT.  When  the  number  of  beams  is  greater  than  7,  the  class-solution  becomes  less  meaningful 
but  may  be  useful  to  “streamline”  clinical  practice.  Finally,  we  mention  that  in  a  sense  any 
optimization  algorithm  has  a  certain  degree  of  subjectivity  due  to  the  difference  in  the  objective 
function  employed  to  model  the  system.  It  would  be  interesting  to  explore  such  model  dependence 
and  further  examine  the  validity  of  the  proposed  class-solution  using  different  model  systems. 
Given  the  fact  that  IMRT  is  being  increasingly  employed  for  prostate  cancer  irradiation,  we 
believe  that  this  type  of  study  may  have  practical  implication  for  clinical  IMRT,  especially  before 
the  users  are  provided  with  a  robust  beam  orientation  optimization  tool  in  the  clinical  IMRT 
treatment  planning  systems. 
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FIGURES 


Dose  (Gy) 


Figure  1  DVH  spectra  of  prostate  (red),  rectum  (green),  and  bladder  for  a  representative  patient  for  all  the 
non-dominated  plans  in  the  Pareto  front. 


Figure  2  Angular  distribution  of  beams  selected  by  the  optimizer,  for  plans  with  3,5-9  beams.  The  length  of  the 
radial  segment  in  the  above  graphs  represents  the  frequency  of  occurrence  of  the  corresponding  direction. 
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Figure  4  Distributions  of  beam  angles  for  3, 5-9  beams,  respectively.  The  colored  lines  represent  directions 
found  in  individual  cases  and  red  bold  lines  represent  the  directions  identified  as  the  class-solutions. 


17 


Figure  5  Comparison  of  plans  with  customized  beams  directions  and  class-solution  directions  a)  Pareto  fronts 
obtained  by  the  two  methods  b)  DVHs  of  class-solution  plan  (broken  line)  and  optimized  plan. 


Figure  6  Comparison  of  class-solution  plans  and  optimised  plans  for  9  beams,  a)  Pareto  fronts  b)  DVHs 
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Figure  7  Color  wash  dose  distributions  in  axial,  coronal,  and  sagittal  planes,  for  the  following  cases:  a)  5  beams 
plan  with  optimized  directions;  b)  5  beams  plan  with  class-solution  directions;  c)  9  beams  plan  with  optimised 
directions;  and  d)  9  beams  plan  with  class-solution  directions.  Structures  are  represented  as  lines  :  rectum 
(green)  bladder  (violet)  and  PTV  (red) 
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TABLES 


Organ 

Max  dose  Dmax 

Critical  dose  DCritic 

Limit  volume  VCritic 

(Gy) 

(Gy) 

(%  of  organ  volume) 

NT 

75 

40 

20 

Left  Femur 

68 

40 

20 

Rectum 

72 

53 

25 

Right  Femur 

68 

40 

20 

Table  I  Dose  limits  and  volume  limits  used  as  constraints  for  OARs.  The  prescribed  dose  was  72  cGy  to  the 
PTV  and  0  cGy  to  the  OARs.  For  each  OAR,  two  constraints  are  imposed:  (a)  the  maximum  dose;  and  (b)  a 
specified  volume  should  not  receive  a  dose  greater  than  the  limit  dose. 
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Purpose;  Clinical  IMRT  treatment  plans  are  currently  done  using  dose-based 
optimization  algorithms,  which  do  not  consider  the  nonlinear  dose-volume  effects  for 
tumors  and  normal  structures.  The  choice  of  structure  specific  importance  factors 
represents  an  additional  degree  of  freedom  of  the  system  and  makes  rigorous 
optimization  intractable.  The  purpose  of  this  work  is  to  circumvent  the  two  problems  by 
developing  a  biologically  more  sensible  yet  clinically  practical  inverse  planning 
framework. 

Methods  and  Materials:  The  dose-volume  status  of  a  structure  was  characterized  by 
using  the  effective  volume  in  the  voxel  domain.  A  new  objective  function  was 
constructed  with  incorporation  of  the  volumetric  information  of  the  system  so  that  the 
figure  of  merit  of  a  given  IMRT  plan  depends  not  only  on  the  dose  deviation  from  the 
desired  dose  distribution  but  also  the  dose-volume  status  of  the  involved  organs.  To 
automate  the  determination  of  the  structure  specific  importance  factors,  we  wrote  the 
conventional  importance  factor  of  an  organ  into  a  product  of  two  components:  (i)  a 
generic  importance  that  parameterizes  the  relative  importance  of  the  organs  in  the  ideal 
situation  when  the  goals  for  all  the  organs  are  met;  (ii)  a  dose-dependent  factor  that 
quantifies  our  level  of  clinical/dosimetric  satisfaction  for  a  given  plan.  The  generic 
importance  can  be  determined  a  priori,  and  in  most  circumstances,  does  not  need 
adjustment,  whereas  the  second  one,  which  is  responsible  for  the  intractable  behavior  of 
the  tradeoff  seen  in  conventional  inverse  planning,  was  determined  automatically.  For 
this  purpose,  we  heuristically  expressed  the  second  factor  as  a  function  of  TCP  or  NTCP. 
After  beam  optimization,  the  TCP/NTCP  was  computed  and  its  value  determined  whether 
the  corresponding  importance  factor  should  be  increased  or  decreased.  The  procedure 
proceeded  in  an  iterative  fashion  until  the  overall  objective  function  stopped  improving. 
An  inverse  planning  module  based  on  the  proposed  formalism  was  implemented  and 
applied  to  a  prostate  case  and  a  head-neck  case. 

Results:  The  incorporation  of  clinical  knowledge  allows  us  to  obtain  IMRT  plans  that 
would  otherwise  be  unattainable  and  makes  it  possible  to  auto-select  the  importance 
factors,  greatly  facilitating  the  inverse  planning  process.  Comparison  of  the  newly 
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proposed  approach  with  the  conventional  inverse  planning  technique  indicated  that,  for 
the  same  target  dose  coverage,  the  critical  structure  sparing  was  substantially  improved 
for  both  cases.  Conversely,  the  dose  to  the  target  volume  can  be  escalated  greatly  while 
maintaining  the  radiation  toxicity  at  the  current  level. 

Conclusions:  The  modeling  of  inverse  planning  objective  function  plays  an  essential  role 
in  the  success  of  IMRT  treatment  and  the  incorporation  of  the  existing  clinical  endpoint 
data  into  inverse  planning  affords  an  effective  way  to  improve  the  sub-optimal 
performance  of  current  IMRT  planning  algorithms.  The  new  formalism  proposed  in  this 
paper  also  reveals  the  relationship  between  different  inverse  planning  schemes  and  sheds 
important  insight  into  the  problem  of  therapeutic  plan  optimization. 


Key  word:  IMRT,  inverse  planning,  dose  optimization,  objective  function,  NTCP 
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Introduction 

An  important  issue  in  inverse  treatment  planning  is  how  to  formalize  the  clinical 
goals  to  objectively  evaluate  the  figures  of  merit  of  different  IMRT  plans.  Despite  of 
intense  research  effort  in  modeling  the  clinical  decision-making  strategies(l-6), 
appropriate  form  of  the  objective  function  for  clinical  IMRT  planning  remains  illusive. 
Presently,  two  types  of  objective  functions  are  being  widely  used:  dose-  or  dose  volume 
histogram  (DVH)-based  (physical  objective  functions)(7-16)  and  dose  response-based 
objective  functions  (biological  objective  functions)(  17-23).  The  underlying  difference 
between  these  models  lies  in  what  endpoints  are  used  to  evaluate  the  treatment  or  which 
fundamental  quantities  are  used  to  define  the  optimality.  The  physical  approach 
emphasizes  the  difference  between  the  calculated  and  prescribed  doses  and  does  not 
consider  the  nonlinear  effects  for  tumors  and  normal  structures.  Dose-volume  constraints 
are  often  introduced  to  select  a  solution  with  certain  shapes  of  the  DVHs  for  the  target 
and  sensitive  structures.  However,  it  is  important  to  note  that  the  construction  of  the  DVH 
constraints  is  a  priori  in  nature.  The  use  of  constraints  can  only  passively  restrict  the 
accessible  DVHs  of  the  final  solution  by  narrowing  the  solution  space  and  the  figures  of 
merits  of  the  physically  realizable  plans  are  not  changed  as  long  as  they  satisfy  the 
constraints.  To  reflect  our  preference  over  certain  DVHs  for  a  structure  it  is  necessary  to 
express  the  objective  of  the  structure  as  a  function  of  the  volumetric  status,  which  has  not 
been  achieved  in  inverse  treatment  planning  up  to  this  point.  On  the  other  hand, 
biological  model-based  optimization  proponents  argue  that  plan  optimization  should  be 
guided  by  estimates  of  biological  effects,  which  depend  on  the  radiation  dose  through  the 
dose  response  function.  The  treatment  objective  in  biological  model-based  inverse 
planning  is  usually  stated  as  the  maximization  of  the  tumor  control  probability  (TCP) 
while  maintaining  the  normal  tissue  complication  probability  (NTCP)  to  within 
acceptable  levels(5,  17-21,  23).  In  principle,  the  biologically  based  models  are  most 
relevant  for  radiotherapy  plan  ranking(24-29).  However,  the  dose-response  function  of 
various  structures  is  not  sufficiently  understood  and,  at  this  point,  there  is  considerable 
controversy  about  the  models  for  computing  dose-response  indices  and  their  use  in 
optimization.  In  reality,  the  use  of  dose-response  indices  for  optimization  may  lead  to 
very  inhomogeneous  target  dose  distributions^,  18-21,  23).  Furthermore,  it  is  difficult 
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for  clinicians  to  specify  the  optimization  criteria  in  terms  of  certain  dose  response  indices 
(e.g.,  TCP,  NTCP  and  P+).  This  issue  becomes  compounded  when  two  or  more 
independently  optimized  plans  are  to  be  combined.  Because  of  these  problems,  the  use  of 
biological  model-based  dose  optimization  has  mainly  been  limited  in  research  community 
and  little  effort  has  been  made  to  implement  the  model  into  commercial  IMRT  planning 
system.  Given  the  fact  that  biological  outcome  is  the  ultimate  endpoint  of  radiation 
therapy,  the  importance  of  the  biological  modeling  and  the  biological  model-based 
inverse  planning  can  never  be  over  stated. 

To  pin  down  the  underlying  problem  of  the  current  inverse  planning  formalism 
and  illustrate  the  need  for  a  clinically  more  relevant  approach,  let  us  take  parotid  glands 
as  an  example.  It  is  well  known  that  the  clinical  end  point  is  the  same  if  the  glands  are 
irradiated  15Gy  to  67%,  or  30Gy  to  30%,  or  45Gy  to  24%  of  the  total  volume(30).  If  a 
dose-based  metric,  such  as  the  commonly  used  quadratic  objective  function,  is  used,  the 
rankings  for  the  three  different  scenarios  would  be  completely  different.  Even  with  the 
use  of  dose-volume  constraints,  it  is  difficult,  if  not  impossible,  to  incorporate  this  type  of 
knowledge  to  correctly  model  the  behavior  of  the  organ  in  response  to  radiation.  Indeed, 
a  constraint  in  optimization  acts  as  a  “boundary  condition”  during  the  optimization  and 
does  not  change  the  rankings  of  dosimetrically  different  plans.  This  example  clearly 
reveals  the  inadequacy  of  the  conventional  dose-based  objective  function  and  suggests 
the  urgent  need  for  a  clinically  more  sensible  model.  Obviously,  a  minimum  requirement 
for  the  model  is  that  it  should  be  consistent  with  the  existing  clinical  outcome  data.  For 
parotid  glands,  for  instance,  the  three  different  DVHs  mentioned  above  should  be  scored 
equally  by  the  objective  function.  We  believe  that  this  type  of  “degeneracy”  is  achieved 
by  effectively  integrating  clinical  end  point  data  into  the  inverse  planning  formulation. 
For  a  given  patient,  as  for  which  of  the  degenerate  DVHs  will  be  selected,  it  will  be 
determined  by  the  optimization  algorithm  with  consideration  of  the  dosimetric/clinical 
requirements  of  other  structures.  The  specifics  of  the  plan  selection  process  will,  of 
course,  depend  on  the  geometric  and  dosimetric  details  of  the  given  patient.  The  example, 
however,  underscores  the  important  role  of  the  existing  clinical  data  in  inverse  planning 
and  emphasizes  the  essential  ingredients  for  a  clinically  realistic  objective  function. 
Towards  developing  a  biologically  more  sensible  yet  clinically  practical  inverse  planning 
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formalism,  in  the  following  we  propose  a  method  to  incorporate  clinical  endpoint  data 
into  the  construction  of  the  IMRT  objective  function  and  attempt  to  bridge  the  gap 
between  the  clinical  decision-making  process  and  the  computational  modeling.  Our  study 
indicates  that  the  clinical  knowledge-based  model  developed  in  this  work  allows  us  to 
more  objectively  rank  radiation  therapy  plans  according  to  their  clinical  merits  and  makes 
it  possible  to  obtain  truly  optimal  IMRT  plans  with  much  reduced  efforts. 


Methods  and  Materials 

Dose-volume  effect 

To  express  the  treatment  objective  of  an  organ  as  a  volumetric  function,  a  critical 
step  is  to  identify  the  quantity  that  characterizes  the  dose-volume  status  of  the  organ  and 
differentiates  one  DVH  state  from  another.  Once  the  quantity  is  selected,  the  objective 
function  can  be  constructed  in  the  representation  defined  by  the  quantity  with 
consideration  of  clinical  outcome  data.  It  appears  that  there  are  two  ways  to  proceed.  One 
is  to  work  in  the  dose-volume  domain,  in  which  the  dose-volume  status  of  an  organ  is 
characterized  by  the  shape  of  the  DVH  curve  itself.  Since,  at  this  point,  clinical  outcome 
data  are  sparse  and  underdetermined,  it  is  impossible  for  us  to  objectively  rank  all 
physically  realizable  plans  based  on  the  data,  necessitating  an  interpolation/extrapolation 
scheme  for  plan  ranking  if  one  chooses  to  work  in  the  dose-volume  domain.  A  more 
convenient  representation  is  perhaps  the  voxel  domain  since  it  is  relatively  intuitive  and 
much  work  has  already  been  done  in  the  past  to  describe  the  dose-volume  effect.  The  key 
here  is  how  to  “extrapolate”  the  clinical  endpoint  observed  at  an  organ  level  to  an 
expression  at  the  individual  voxel  level.  However,  it  is  important  to  note  that  the  dose- 
volume  presentation  is,  in  principle,  a  more  rigorous  domain  to  deal  with  the  dose- 
volume  effect.  The  equivalence  mapped  to  the  voxel  domain  is  generally  established 
phenomenologically  based  on  sparse  clinical  data  and  represents,  at  best,  an 
approximation  of  the  true  dose-volume  response  relation. 

In  this  research  we  work  in  the  voxel  domain.  Generally,  the  dose  response  of  a 
structure  with  respect  to  the  irradiated  dose  and  volume  is  complicated.  The  fact  that  the 
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dose  distribution  in  tumor  or  a  sensitive  structure  is  generally  inhomogeneous  makes  the 
establishment  of  such  a  relationship  even  more  intractable.  Over  the  last  two  decades, 
attempts  have  been  made  by  many  researchers  to  capture  the  main  feature(s)  of  the  dose 
volume  effects.  A  power  law  model  represents  one  of  the  successful  techniques  in  dealing 
with  the  dose-volume  effects  of  sensitive  structures(28,  31).  In  this  model  an  equivalent 
dose  uniformly  irradiating  the  whole  organ,  Deq,  can  be  used  to  represent  the  situation  in 
which  a  fractional  partial  volume,  v,  is  irradiated  to  a  dose,  D,  by  a  simple  power  law 
model:  Deq=v1/nD.  A  remarkable  characteristic  of  this  model  is  that,  although  only  a 
single  organ-specific  parameter,  n,  is  used,  clinical  and  biological  data  has  shown  that 
this  power  law  holds  well  at  low  complication  levels(28,  32).  Based  on  this  relation, 
Mohan  et  al(17)  introduced  the  concept  of  effective  dose  to  represent  a  non-uniform  dose 
distribution  in  a  sensitive  structure.  Kutcher  and  Burman(29)  applied  the  same  power 
model  independently  to  each  volume  element  of  the  histogram  and  introduced  the 
concept  of  effective  volume  to  reduce  the  DVH  of  an  inhomogeneous  dose  distribution  in 
a  sensitive  structure  to  a  uniform  dose  distribution.  Following  their  study,  we  define  the 
effective  volume  (A  Veff)i  for  a  voxel  i  with  volume  A  F  and  dose  Dc(i)  as  follows 
(AFejr)i=AF(Dl(i)/Dref)Un  (1) 

and  extend  this  concept  to  handle  the  voxels  in  the  tumor  target,  where  n  is  an  organ- 
dependent  parameter  and  Dref  is  the  reference  dose.  For  a  sensitive  structure,  n  is  a  small 
positive  number  (0<n<l)  and  the  value  of  parameter  n  reflects  the  architecture  (serial  or 
parallel)  of  the  sensitive  structure.  For  a  target,  n  should  be  assigned  with  a  small 
negative  value  (-l<n<0).  The  biological  meaning  of  Eq.  (1)  is  that  for  a  sensitive 
structure  a  small  volume  receiving  a  higher  dose  than  a  reference  dose  would  be 
equivalent  to  a  larger  volume  receiving  the  reference  dose;  for  a  target,  a  volume  with  a 
lower  dose  would  have  a  larger  effective  volume.  The  effective  volumes  of  all  voxels 
reflect  the  DVH  status  of  the  given  organ,  and  for  inverse  planning,  this  permits  us  to 
deal  with  the  complicated  dose-volume  effect  in  the  voxel  domain. 

Dose-volume  based  objective  function 

The  objective  function  expressed  as  a  function  of  the  effective  volume  in  the 
voxel  domain  should  take  the  form  of 
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(2) 


L=fM  a*V),}), 

Generally,  the  dose-volume  effect  suggests  that  the  voxels  receiving  different  doses  are 
inequivalent:  the  one  with  a  larger  effective  volume  (higher  dose  for  a  critical  organ) 
should  be  penalized  more  when  compared  to  a  voxel  with  a  smaller  effective  volume 
(lower  dose).  Thus  we  heuristically  write  the  fa  in  the  following  form 

-|2 

/*=i+^1Zr/(MVX  +,72  +•••’  (3) 

i  L  I 

where  r,  is  the  importance  factor  of  the  z'-th  voxel,  representing  the  intrastructural  tradeoff 
due  to  physical/clinical  requirements  other  than  the  dose-volume  based  penalty,  r\i  and  r\2 
are  phenomenological  parameters  of  the  model.  In  Eq.  (3),  the  third  (and  higher  order) 
term  emphasizes  more  on  the  voxels  with  high  effective  volumes,  whereas  the  first  and 
second  terms  ensure  that  the  voxels  with  low  effective  volumes  receive  an  adequate 
penalty.  We  typically  set  rt=  1  in  Eq.  (3),  unless  there  are  other  physical/clinical 

considerations  (e.g.,  when  the  density  of  clonogenes  varies  spatially  (33)).  In  this  work, 
unless  specifically  mentioned,  we  set  tji  =land  z/2=z/3=...  =0. 


Hybrid  of  dose-based  and  dose  volume-based  objective  functions 

Equation  (3)  provides  a  good  description  of  the  dose-volume  effect.  With  proper 
choice  of  the  parameter,  n,  the  clinically  observed  dose  volume  effect  can  be  reproduced 
by  the  objective  function.  In  reality,  other  clinical  requirements,  such  as  the  target  dose 
homogeneity,  should  also  be  considered.  A  more  general  form  of  inverse  planning 
objective  function  can  be  written  as  a  hybrid  of  the  dose-volume  based  and  the  dose- 
based  functions.  In  this  situation,  the  overall  objective  function  of  the  system  takes  the 
form  of 


F  =  t,r—^{\  +  n\[Dc(j)IDTtKfyn' 


}\Dc(i)-DT0(i)\' 


+ 2/ ^  tt  2  A  ■ +  <  [A  (0  /  }  Dc(i)k* , 


where  Eq.  (1)  has  been  incorporated,  rt  and  ra  are  the  structure  specific  importance  factor 
for  target  x  or  sensitive  structure  o,  tx  and  sa  are  the  numbers  of  targets  and  sensitive 
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structures,  nz  and  n„  are  the  n  parameters  for  target  x  and  sensitive  structure  <r,  Dx,ref  and 
Da, ref  are  the  reference  doses  for  target  x  and  sensitive  structure  <7,  respectively.  The  factor 

|/)c(/)-£)0r(/)| '  for  target  or  Dc(i)k°  for  a  sensitive  structure  represents  the  contribution 

from  dosimetric  deviation  from  the  ideal  situation.  If  the  kr  and  kB  are  set  to  zero,  Eq.  (4) 
becomes  identical  to  Eq.  (3)  and  the  objective  function  becomes  purely  dose-volume 
driven.  In  particular,  if  we  set  ka  to  zero  and  ^  to  a  non-zero  value,  the  objective  function 
for  a  target  becomes  a  hybrid  of  dose-volume  and  dose  based,  whereas  the  objective 
functions  for  critical  structures  remain  to  be  purely  dose-volume  based.  On  the  other 
hand,  when  all  the  n  parameters  in  Eq.  (4)  are  set  to  be  +oo,  no  dose-volume  effects  are 
considered  and  Eq.  (4)  is  reduced  to  the  conventional  dose-based  objective  function.  It  is 
interesting  to  note  that,  according  to  Eqs.  (2)  and  (4),  the  added  dose-based  factor  can  be 
interpreted  as  a  voxel-dependent  importance  factor  which  increases  with  the  dosimetric 
deviation  from  the  ideal  situation  for  the  target  or  a  sensitive  structure.  It  seems  to  make 
intuitive  sense  to  assign  a  high  weight  (or  importance)  to  those  voxels  that  receive  doses 
farther  away  from  the  ideal  situations  so  that  more  penalty  can  be  applied  to  these  voxels. 
Conversely,  the  dose-volume  factor,  fa,  can  be  regarded  as  the  voxel-dependent 

importance  factor  for  a  dose-based  objective  function.  The  local  importance  factor  now 
increases  with  the  effective  volume  at  the  voxel.  We  emphasize  that  these  voxel-based 
importance  factors  are  derived  directly  from  the  existing  clinical  knowledge  instead  of 
heuristic  arguments. 

Automatic  determination  of  structure  specific  importance  factors 

The  selection  of  structure  specific  importance  factors,  rx  or  ra  in  Eq.  (4),  is 
generally  done  empirically  with  trial-and-error.  Here  we  describe  an  alternative  and  more 
automated  approach  for  solving  the  problem  and  show  that  the  selection  of  structure 
specific  importance  factors  can  be  facilitated  by  using  computer  optimization.  The  key  of 
success  is  to  establish  an  effective  method  to  express  the  structural  importance  factor  in 
terms  of  physically  or  clinically  more  meaningful  quantities.  For  this  purpose  we  write 
the  importance  factor  of  an  involved  structure  into  a  product  of  two  components:  (i)  a 
generic  importance  that  parameterizes  the  relative  importance  of  the  organs  in  an  ideal 
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situation  when  the  goals  for  the  organs  are  met;  and  (ii)  a  dose-dependent  factor  that 
quantifies  our  level  of  clinical/dosimetric  satisfaction  for  a  given  plan.  The  first  factor  can 
be  determined  a  priori,  and  in  most  circumstances,  does  not  need  adjustment  (generally 
speaking,  the  value  of  r/  is  determined  based  on  the  treatment  modality  and  the  patient’s 

overall  condition,  age,  types  of  complications  expected,  and  so  on),  whereas  the  second 
one  is  responsible  for  the  intractable  behavior  of  the  tradeoff  in  conventional  planning 
and  can  be  automatically  determined.  This  decomposition  is  essentially  to  normalize  the 
conventional  importance  factor  in  terms  of  our  clinical  goal  for  the  structure  under 
discussion.  Because  of  this  decomposition,  the  meaning  of  the  importance  factor  is  more 
transparent  and  the  determination  of  the  factors  becomes  straightforward. 
Mathematically,  we  write  the  organ  specific  importance  factor  of  a  sensitive  structure  as 
ra  =  r/r* ,  where  r/ represents  the  first  contribution  described  above  (the  desired 

weighting  among  different  structures  in  an  ideal  situation),  and  r*  is  the  second 
component  and  is  defined  as  a  relative  ratio  against  a  known  clinical  end  point  (for 
example,  5%  complication  rate  for  a  sensitive  structure).  In  this  study  r/  was  set 
empirically  (See  tables  II  and  III  for  examples),  r^  is  updated  according  to  the  DVH  or 
the  dose  distribution  during  the  optimization  process  and  reflects  the  most  current  status 
of  tradeoff  in  the  system.  For  a  sensitive  structure,  the  value  of  r*  is  defined  as,  in  this 

study,  the  ratio  of  NTCP  against  a  desired  complication  probability,  say  5%.  Generally, 
the  importance  of  a  sensitive  structure  should  be  increased  in  next  iteration  if  NTCP  is 
high,  and  vice  versa.  We  found  that  a  simple  linear  relation  between  r*  and  NTCP, 

r%  =  NTCPa  +  8  (5) 

describes  the  tradeoff  behavior  of  the  system  well,  where  NTCPa  is  the  NTCP 
corresponding  to  the  dose  distribution  at  the  current  iteration  for  structure  a,  d  is  a  cutoff 
factor  for  NTCP,  which  is  introduced  to  ensure  the  sensitive  structure  receiving  a 
minimum  penalty  even  if  its  NTCP  is  close  to  zero.  We  set  d  as  an  organ-independent 
constant  of  0.01%. 

NTCP  was  assessed  using  Lyman’s  model  in  this  study.  For  non-uniform 
irradiation,  the  Kutcher-Burman  effective-volume  DVH  reduction  method(29)  is  used  to 
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transform  a  DVH  into  a  uniform  irradiation  on  an  effective  partial  volume.  There  are 
three  parameters  in  this  model  (n,  m,  TD50/5).  The  parameter  n  represents  the  volume 
effect  and  the  parameter  m  is  the  steepness  of  the  dose  complication  curve  for  a  fixed 
partial  volume.  Model  parameters  used  in  this  study  were  those  fitted  to  the  model  by 
Burman  et  al(34)  for  the  normal  tissue  tolerance  data  compiled  by  Emami  et  al(35). 


Computational  Algorithm 

After  considering  the  automatic  tradeoff  strategy,  the  generalized  objective 
function  takes  the  following  form 

r=l  ^  ’  t  1=1 


We  implemented  a  software  module  to  optimize  the  objective  function  (6)  in  the  platform 
of  the  PLUNC  treatment  planning  system  (University  of  North  Carolina,  Chapel  Hill, 
NC).  The  dose  calculation  engine  and  varieties  of  evaluation  tools  of  the  PLUNC  system 
were  used  to  evaluate  and  compare  the  optimization  results.  The  calculation  was 
performed  on  a  PC  with  P4  1.7GHz  and  1024MB  RAM.  The  ray-by-ray  iterative 
algorithm  (SIITP)  reported  earlier(9,  36)  was  employed  to  obtain  the  optimal  beam 
intensity  profiles.  The  values  of  kx  and  ka  in  Eq.  (6)  were  set  to  be  2,  but  the  behavior  of 
the  system  for  a  few  other  combinations  of  K  and  ka  were  also  checked  for  the  prostate 
case.  The  reference  dose,  Da<ref,  was  chosen  to  be  TD5/5  of  the  corresponding  critical 
organ.  For  the  target,  DT,re/was  set  as  the  prescription  dose.  Figure  1  shows  the  flow  chart 
of  the  calculation  process.  In  current  study  we  specify  a  maximum  number  of  iterations 
as  the  termination  condition  of  the  optimization  process.  The  DVHs  can  be  inspected  in 
each  iterative  step  to  visually  monitor  the  optimization  process. 


Case  studies 

Two  cases,  a  prostate  case  and  a  head-neck  case,  were  used  to  evaluate  the 
proposed  inverse  planning  formalism.  The  optimization  results  were  compared  with  those 
obtained  using  the  conventional  dose-based  optimization  method,  which  was  described  in 
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details  by  Xing  et  al  (9,  36).  The  optimization  parameters  in  the  dose-based  method  were 
adjusted  by  trial-and-error  to  obtain  an  “optimal”  plan. 

In  the  prostate  IMRT  case,  the  target  volume  included  the  prostate  and  seminal 
vesicles.  The  sensitive  structures  involved  in  this  study  were  rectum,  bladder  and  femoral 
heads.  All  the  IMRT  plans  used  identical  configuration  of  five  equally  spaced  15MV 
photon  beams  with  gantry  angles  of  0°,  72°,  144°,  216°,  and  288°  (in  IEC  convention). 
The  plans  were  normalized  to  deliver  the  prescription  dose  of  70Gy  to  99%  of  the  target 
volume.  The  parameter  nT  was  chosen  to  be  -0.2  for  the  target.  The  parameters  used  for 
the  computation  of  the  NTCPs  for  rectum,  bladder  and  femoral  heads  are  listed  in  Table 
I,  which  were  obtained  by  Emami  at  al(35)  and  Burman  et  al(34).  Table  II  summarizes 
the  optimization  parameters  for  both  the  newly  proposed  and  dose-based  approaches. 

For  the  prostate  case,  we  also  studied  the  influence  of  two  more  combinations  of 
ka  and  kT.  These  included  {ka= 2,  kT= 4)  and  (ka= 0,  kT=2).  In  the  latter  case,  we  have 
included  a  higher  order  term  of  the  dose-volume  effect  (the  third  term  in  Eq.  (3)  with 
t|2=1)  to  ensure  that  the  high  effective  volume  voxels  are  penalized  enough  in  the  absence 
of  the  dose-based  factor. 

In  the  head-and-neck  case,  the  organs  at  risk  included  the  eyes,  optic  nerves,  optic 
chiasm,  brainstem,  spinal  cord  and  parotids.  Two  treatment  targets  were  included  in  this 
case:  the  gross  target  volume  (GTV)  and  the  clinical  target  volume  (CTV),  which 
includes  the  microscopic  disease  region  surrounding  the  GTV.  The  plan  was  normalized 
to  deliver  a  prescription  dose  of  70Gy  to  at  least  99%  of  the  GTV  and  62Gy  to  at  least 
95%  of  the  CTV.  Nine  equally  spaced  6MV  coplanar  beams  (0°,  40°,  80°,  120°,  160°, 
200°,  240°,  280°,  and  320°  in  IEC  convention)  were  used  for  this  case.  The  parameter  nT 
were  -0.5  for  both  GTV  and  CTV.  The  parameters  used  for  the  computation  of  the 
NTCPs  of  the  sensitive  structures  are  also  obtained  from  the  same  source  stated  earlier 
and  are  listed  in  Table  I.  The  optimization  parameters  for  both  the  two  techniques  are 
summarized  in  Table  III. 


Results 

Prostate  IMRT  plans 
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Figures  2  and  3  summarize  the  results  of  the  two  IMRT  plans  obtained  using  the 
newly  proposed  and  conventional  techniques.  Figure  2  compares  the  isodose  distributions 
in  two  transverse  slices  and  a  sagittal  slice  for  the  two  plans.  The  DVHs  of  the  prostate 
target  and  sensitive  structures  are  plotted  in  figure  3,  in  which  the  solid  and  dashed  lines 
represent  the  DVHs  obtained  using  the  new  and  conventional  approaches,  respectively. 
The  calculated  NTCPs  of  rectum,  bladder  and  femoral  heads  for  both  IMRT  plans  are 
listed  in  table  IV.  According  to  the  table,  it  is  seen  that  the  NTCPs  of  the  sensitive 
structures  are  improved  significantly.  For  the  rectum,  for  example,  the  NTCP  is  reduced 
from  0.45%  to  0.03%.  Our  results  also  indicate  that  the  main  compromise  in  a  prostate 
IMRT  treatment  seems  to  be  between  the  tumor  coverage  and  the  rectum  complication 
because  the  NTCP  of  the  rectum  is  much  higher  than  that  of  other  sensitive  structures. 

The  above  results  demonstrate  that,  for  comparable  target  coverage,  the  new 
inverse  planning  technique  greatly  improves  the  critical  structure  sparing,  especially  the 
rectum  sparing.  By  comparing  the  isodose  distributions  of  the  two  plans  (figure  1),  it  is 
seen  that  the  dose  gradient  at  the  interface  between  the  target  and  the  rectum  is  much 
steeper  for  the  IMRT  plan  obtained  with  the  new  formalism.  Furthermore,  it  is  intriguing 
that  the  non-sensitive  structure  normal  tissue  also  receives  less  dose  in  comparison  with 
that  of  the  dose-based  optimization.  Our  results  suggest  that  the  improvement  in  the 
critical  structure  sparing  is  achieved  not  at  the  cost  of  higher  target  dose  inhomogeneity, 
which  is  commonly  seen  in  IMRT  plan  optimization. 

The  resultant  DVHs  when  ka=2  and  kT  was  increased  from  2  to  4  in  Eq.  (6)  are 
plotted  in  Fig.  3  as  the  dotted  curves.  While  the  target  dose  uniformity  is  improved  when 
the  kT  increases,  the  doses  to  the  rectum  and  bladder  are  worsened.  The  results  make 
intuitive  sense  as  when  the  kT  increases,  more  penalty  is  applied  toward  dosimetric 
deviation  from  the  prescription.  The  DVHs  when  the  objective  function  for  the  target  is  a 
hybrid  of  dose-volume  and  dose  based  functions  (kT=  2)  and  that  for  the  sensitive 
structures  are  purely  dose-volume  based  (ka= 0)  are  shown  in  Fig.  3  as  dash-dotted  curves. 
In  this  case,  a  high  order  term  of  the  dose-volume  effect  in  Eq.  (3)  was  added  to  ensure 
that  the  high  effective  volume  voxels  are  penalized  enough  in  the  absence  of  the  dose- 
based  factor.  Interestingly,  as  can  be  seen  from  Fig.  3,  the  results  so  obtained  were  very 
similar  to  that  obtained  with  the  hybrid  objective  function. 
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Head-and-neck  IMRT  plans 

The  two  IMRT  plans  obtained  using  the  newly  proposed  and  the  dose-based 
techniques  are  summarize  in  figures  4  and  5.  Figure  4  shows  the  isodose  distributions  in 
three  transverse  slices,  one  sagittal  slice  and  one  coronal  slice  for  the  two  plans  and  figure 
5  compares  the  DVHs  of  the  targets  and  sensitive  structures.  In  figure  5  the  solid  and 
dashed  lines  represent  the  DVHs  obtained  using  the  newly  proposed  and  conventional 
approaches,  respectively.  The  calculated  NTCPs  of  eyes,  optical  nerves,  optical  chiasm, 
brainstem,  spinal  cord  and  parotids  for  both  plans  are  shown  in  table  V.  As  seen  from  the 
isodose  distributions  (figure  3)  and  DVHs  (figure  4),  with  comparable  GTV  and  CTV 
dose  coverage  and  dose  homogeneity,  the  doses  to  the  sensitive  structures  are 
dramatically  reduced.  The  dose  reduction  is  particularly  pronounced  in  the  spinal  cord, 
brainstem,  parotids,  and  eyes.  For  the  left  and  right  parotids,  for  example,  the  fractional 
volume  receiving  a  dose  above  25  Gy  is  reduced  from  35%  to  20%  and  15%, 
respectively.  Consistent  with  the  enhanced  dosimetric  conformality  and  similar  to  the 
prostate  case,  much  steeper  dose  gradient  occurs  near  the  boundary  of  the  target  volume. 
The  dose  to  the  non-sensitive  structure  normal  tissue  is  also  less  in  comparison  with  the 
conventional  IMRT  plan.  While  the  NTCPs  of  the  optical  nerves  and  optical  chiasm  are 
small  and  difficult  to  draw  conclusion,  from  Table  V,  it  is  quite  clear  that  the  NTCPs  of 
the  eyes,  parotids,  spinal  cord  and  brainstem  are  improved  significantly.  We  emphasize 
once  again  that  the  significant  improvement  in  sensitive  structure  sparing  is  achieved 
without  deteriorating  the  dose  coverage  of  the  GTV  and  CTV. 


Discussion 

The  currently  available  dose-based  objective  functions  do  not  truly  reflect  the 
nonlinear  relationship  between  the  dose  and  the  response  of  tumors  and  tissues  and  it  is 
highly  desirable  to  incorporate  clinical  outcome  data  in  the  formulation  of  inverse 
planning  to  guide  the  plan  optimization  process.  While  the  dose  dependence  of  a  clinical 
endpoint  may  be  degenerate  in  the  sense  that  it  may  be  caused  by  a  variety  of  dose 
distributions  or  DVHs,  there  exists  no  mechanism  in  conventional  inverse  planning  to 
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model  the  phenomenon.  The  irradiation  of  parotid  glands  mentioned  in  the  Introduction 
represents  an  example  of  this.  On  the  other  hand,  the  conventional  objective  function  may 
impose  some  unrealistic  degeneracy  that  is  inconsistent  with  clinical  experience.  For 
example,  assume  that  in  a  treatment  plan  the  prescription  dose  to  the  target  is  70Gy  and 
that  the  tumor  is  divided  into  two  parts  with  the  same  volume,  one  receiving  a  dose  60Gy 
and  another  80Gy.  The  penalty  values  of  the  two  scenarios  would  be  the  same  according 
to  the  conventional  quadratic  function.  Obviously,  the  two  different  dose  distributions 
would  lead  to  different  outcome.  The  cold  part,  which  would  greatly  diminish  the  tumor 
control,  is  more  detrimental  than  the  hot  spot. 

In  this  paper  we  have  established  a  general  inverse  planning  framework  in  which 
the  penalty  at  a  voxel  depends  not  only  on  the  dose  deviation  from  the  desired  value  but 
also  the  dose-volume  status  of  the  involved  organs.  The  technique  circumvents  the 
problems  mentioned  above  and  makes  it  possible  to  take  advantage  of  the  clinical  data. 
Our  study  shows  that  the  incorporation  of  existing  clinical  knowledge  can  greatly 
facilitate  the  inverse  planning  process  and  allows  us  to  obtain  IMRT  plans  that  would 
otherwise  be  unattainable.  For  the  same  target  dose  coverage,  we  found  that  the  critical 
structure  sparing  was  substantially  improved  for  both  cases.  Conversely,  the  dose  to  the 
target  volume  can  be  escalated  greatly  while  maintaining  the  radiation  toxicity  at  the 
current  level.  It  is  remarkable  that  such  the  significant  improvement  is  resulted  purely 
from  a  better  modeling  of  the  system.  We  attribute  the  significant  improvement  in  dose 
conformality  to  the  more  adequate  modeling  of  the  system  and  the  enlarged  universe  of 
solution  space  when  the  dose-volume  effects  are  taken  into  account.  Physically,  we 
believe  that  the  superior  performance  of  the  new  formalism  arises  from  the  adequate 
modulation  of  the  voxel  dependent  weighting  induced  by  the  dose-volume  factor /„.  (see 
Eqs.  (3)  and  (4)).  In  conventional  dose-based  objective  function,  fa=  1 ,  and  a  tacit 
assumption  that  all  points  within  a  structure  are  equivalent  has  been  made.  The  use  of 
dose-volume  factor  fa  given  by  Eq.  (3)  enables  us  to  weight  different  voxels  according  to 
the  local  doses.  In  this  way,  we  can  effectively  "boost"  those  target  regions  where  the 
doses  are  low  or  penalize  more  to  those  sensitive  structure  regions  where  the  doses  are 
high.  The  dose-volume  induced  voxel  inequality  is  an  important  feature  of  the  new 
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inverse  planning  formalism  and  is  the  main  driving-force  in  improving  the  dose 
distributions. 

The  use  of  clinical  knowledge  can  also  facilitate  the  determination  of  the  structure 
specific  importance  factors.  While  the  general  influence  of  the  importance  factors  on  the 
solution  is  known,  the  specific  response  of  the  plan  to  a  variation  in  the  factors  is  not 
clear  until  the  dose  optimization  is  done,  which  necessitates  a  manual  trial-and-error 
adjustment  of  the  factors  to  achieve  an  acceptable  tradeoff.  The  underlying  deficiency  of 
the  conventional  approach  is  that  the  importance  factors  are  purely  heuristic  and  lack  of 
physical/clinical  meanings.  In  this  work  we  proposed  a  new  scheme  for  modeling  the 
tradeoff  and  develop  an  algorithm  to  auto-determine  the  factors.  We  wrote  the 
importance  factor  of  an  organ  into  a  product  of  a  generic  importance  and  a  dose- 
dependent  factor.  The  latter  was  related  to  TCP  or  NTCP  according  to  Eq.  (5).  After 
beam  optimization,  the  dose-dependent  factors  were  increased  or  decreased  according  to 
the  values  of  TCPs/NTCPs.  This  procedure  is  similar  to  that  reported  by  Xing  et  al  (37), 
where  a  DVH-based  “distance”  was  used  for  the  assessment  of  the  tradeoff  status  after 
each  optimization.  In  reality,  other  types  of  plan  evaluation  indices,  such  as 
mean/maximum/  minimum  doses,  can  also  be  employed  for  the  purpose.  We  noticed  that, 
with  the  use  of  new  objective  function,  the  final  solution  becomes  much  less  sensitive  to 
choice  of  { #•/ } .  This  characteristic  of  the  new  formalism  may  have  practical  implications 
in  simplifying  the  inverse  planning  process. 

We  should  acknowledge  that  the  new  technique,  just  like  any  other  dose  response- 
based  technique,  may  be  limited  by  the  scarcity  and  uncertainty  of  biological  data  and  the 
limited  predictive  power  of  the  TCP/NTCP  models.  It  is  important  to  note,  however,  that 
the  TCP/NTCP  is  used  as  a  relative  ranking  in  our  plan  optimization  algorithm  instead  of 
a  clinical  decision-making  tool.  Because  of  the  phenomenological  nature  of  the  modeling, 
one  may  further  modify  the  structure  specific  importance  factors  manually  to  achieve  a 
certain  clinical  goal.  The  proposed  technique  can,  at  least,  provide  us  with  a  good  starting 
point  for  the  fine-tuning.  Our  experience,  along  with  the  results  shown  in  the  above 
section,  indicates  that  the  technique  proposed  in  this  work  is  capable  of  generating 
clinically  sensible  plans  and  is  much  more  efficient  than  the  manual  selection  process. 
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Finally,  we  mention  that  the  optimality  of  the  dose  distributions  can  also  be  cast 
into  the  realm  of  equivalent  uniform  dose  (EUD)(20,  21,  27),  which  is  defined  as  the 
biologically  equivalent  dose  which,  if  given  uniformly,  leads  to  the  same  cell  killing  as 
the  actual  non-uniform  dose  distribution.  It  is  interesting  to  note  that  the  EUD-based 
objective  function  can  be  regarded  as  a  special  case  of  the  general  objective  function  (6). 

( i  N  Y/a 

Indeed,  assuming  EUD  =  —  D“  (/)  ,  a  =l/n,  we  can  rewrite  Eq.  (3)  into 

i  J 

F  =  Xrr[l  +  77i  (EUD/ EUDrref)a'  +71'2(EUD/EUDtref)2a’  +...]  + 

r=l 

f>Jl  +  t1s1(EUD/EUD(T^  +  ij^EUD  /  EUD^  +...],  (7) 

<7=1 

which  becomes  a  function  of  EUD.  Different  from  the  EUD-based  model,  the  general 
hybrid  objective  function  given  in  Eq.  (6)  treats  the  dose-volume  effect  at  a  more 
fundamental  voxel  level  with  the  actual  radiation  dose  considered,  which  is  more  flexible 
than  the  EUD  defined  at  a  structure  level.  Because  of  this,  other  clinical/dosimetric 
requirements  can  be  easily  integrated.  Our  study  for  the  prostate  case  suggests  that  it  is 
necessary  to  include  the  higher  order  contribution(s)  if  Eq.  (3)  or  (7)  is  used  to 
appropriately  model  a  sensitive  structure.  Alternatively,  a  hybrid  of  dose-volume  and 
dose  based  objective  function,  as  given  by  Eq.  (6),  can  yield  equally  good  plans.  In 
practice,  equation  (6)  is  quite  broad  and  seems  to  model  the  inverse  planning  system 
effectively.  It  may  also  find  a  natural  application  in  functional  image-guided  IMRT, 
where  the  goal  is  generally  to  produce  a  spatially  inhomogeneous  dose  distribution  (33). 
Finally,  we  note  that  the  formalism  does  not  involve  the  prescription  of  EUD,  which 
could  be  problematic  for  practical  implementation  of  an  EUD-based  model. 


Conclusion 

Inverse  planning  is  an  important  step  in  IMRT  and  its  performance  crucially 
determines  the  quality  of  IMRT  treatment  plans.  In  this  work,  we  provide  a  mechanism 
for  incorporating  clinical  end  point  data  into  inverse  treatment  planning  process  and 
established  a  clinically  practicable  inverse  planning  framework.  We  employed  the 
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effective  volume  to  take  the  volumetric  effects  of  the  involved  organs  into  account.  The 
new  formalism  sheds  important  insight  into  the  problem  of  therapeutic  plan  optimization. 
An  algorithm  for  using  computer  to  aid  the  determination  of  structure  specific  importance 
factors  was  also  developed.  A  key  step  for  accomplishing  the  auto-determination  of  the 
importance  factors  is  the  decomposition  of  the  conventional  importance  factor  into  a 
generic  importance  and  a  dose-dependent  component.  Two  case  studies  were  presented  to 
demonstrate  the  advantages  of  the  proposed  objective  function.  Comparison  of  the  newly 
proposed  approach  with  the  conventional  inverse  planning  technique  indicated  that  the 
algorithm  is  capable  of  greatly  improving  the  sensitive  structure  sparing  with  comparable 
target  dose  coverage  and  homogeneity.  Conversely,  the  dose  to  the  target  volume  can  be 
escalated  greatly  while  maintaining  the  radiation  toxicity  at  the  current  level.  Considering 
the  level  of  improvement  of  the  proposed  technique  and  the  fact  that  it  is  resulted  purely  from 
more  intelligent  computing  algorithm,  as  apposed  to  any  other  exogenous  means  (e.g.,  the 
use  of  radiosensitizer,  whose  effects  to  the  patients  may  not  be  completely  clear),  we  believe 
that  the  potential  impact  of  the  work  on  cancer  treatment  is  significant. 
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Table  I  The  radiological  parameters  for  various  sensitive  structures  used  in  this  study 


Sensitive  structures 

n 

m 

Dso/s  (Gy) 

D5/s  (Gy) 

Bladder 

0.50 

0.11 

80 

65 

Rectum 

0.12 

0.15 

80 

60 

Femoral  head 

0.12 

65 

52 

Eye  lens 

0.27 

18 

10 

Optic  nerve 

0.14 

65 

50 

Optic  chiasm 

0.14 

65 

50 

Spinal  cord 

0.175 

66.5 

47 

Brainstem 

0.14 

65 

50 

0.18 

46 

32 

Table  II  Summary  of  the  optimization  parameters  used  in  the  dose-based  and  proposed 

approaches  for  the  prostate  case 


The  dose-based  approach 

The  proposed 
approach 

Organs 

Relative 

importance 

factors 

Target  prescription 
and  OAR  tolerance 
doses  (Gy) 

Generic  importance 
factors  (r/) 

Target 

5.0 

78 

5 

Bladder 

1.2 

48 

2 

Rectum 

1.8 

43 

2 

Femoral  head  (R) 

1.0 

32 

1 

Femoral  head  (L) 

1.0 

32 

1 

Normal  tissue 

0.5 

65 

0.3 
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Table  III  Summary  of  the  optimization  parameters  used  in  the  dose-based  and  proposed 

approaches  for  the  Head-and-neck  case 


The  dose-based  approach 

The  proposed 
approach 

Organs 

Relative 

importance 

factors 

Target  prescription 
and  OAR  tolerance 
doses  (Gy) 

Generic  importance 
factors  (r/) 

GTV 

3.0 

70 

4.0 

CTV 

4.0 

62 

6.0 

Spinal  cord 

2.0 

30 

3.0 

Brainstem 

1.5 

30 

2.0 

Left  optic  nerve 

1.0 

25 

1.0 

Right  optic  nerve 

1.0 

25 

1.0 

Left  Eye 

2.0 

6 

3.0 

Right  Eye 

2.0 

6 

3.0 

Left  parotid 

1.2 

25 

1.0 

Right  parotid 

1.2 

25 

1.0 

Optic  chiasm 

1.0 

25 

1.0 

Normal  tissue 

0.5 

40 

0.5 

Table  IV  Comparison  of  the  normal  tissue  complication  probabilities  (NTCP)  for  the  two 

IMRT  plans  for  the  prostate  case 


NTCP  (%) 

The  dose-based  IMRT  plan 

The  proposed  IMRT  plan 

Bladder 

0.017 

0.00030 

Rectum 

0.45 

0.029 

Femoral  head  (R) 

0.000076 

0.0000038 

Femoral  head  (L) 

0.000032 

0.000015 
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Table  V  Comparison  of  the  normal  tissue  complication  probabilities  (NTCP)  for  the  two 

IMRT  plans  for  the  head-and-neck  case 


NTCP (%) 

The  dose-based  IMRT  plan 

The  proposed  IMRT  plan 

Spinal  cord 

0.043 

0.0025 

Brainstem 

0.012 

0.0040 

Left  Eye 

0.27 

0.18 

Right  Eye 

0.24 

0.12 

Left  parotid 

0.21 

0.056 

Right  parotid 

0.22 

0.064 

Optic  chiasm 

0.00024 

0.00064 

Left  optic  nerve 

0.000064 

0.0000075 

Right  optic  nerve 

0.000043 

0.000025 

25 


Legends 


Fig.  1.  A  flow  chart  of  the  proposed  optimization  process. 

Fig.  2.  Comparison  of  the  isodose  distributions  of  the  two  prostate  IMRT  plans:  (a)  the 
conventional  dose-based  approach;  (b)  the  newly  proposed  approach.  The  results  on  two 
transverse  slices  and  a  sagittal  slice  are  shown. 

Fig.  3.  Comparison  of  Dose  Volume  Histograms  (DVHs)  of  the  prostate  IMRT  plans 
obtained  using  the  proposed  approach  (solid  curves)  and  the  conventional  dose-based 
approach  (dash  lines).  The  dotted  curves  represent  the  results  obtained  with  £,=4  and 
kff=2  in  Eq.  (6).  The  dash-dotted  curves  are  the  DVHs  with  1q=2  and  k„= 0  (a  higher  order 
term,  the  third  term  in  Eq.  (3),  was  included  during  the  optimization). 

Fig.  4.  Comparison  of  the  isodose  distributions  of  the  two  head-and-neck  IMRT  plans:  (a) 
the  conventional  dose-based  approach;  (b)  the  newly  proposed  approach.  The  results  on 
three  transverse  slices,  one  sagittal  slice  and  one  coronal  slice  are  shown. 

Fig.  5.  Comparison  of  the  Dose  Volume  Histograms  (DVHs)  of  the  two  head-and-neck 
IMRT  plans  obtained  using  our  newly  proposed  approach  (solid  curves)  and  the 
conventional  dose-based  approach  (dash  lines). 
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Abstract 


The  detection  of  Citrate  at  main  magnetic  field  of  3  Tesla  is  investigated.  Citrate 
is  an  important  metabolite  often  used  to  aid  the  detection  of  Prostate  Cancer  but 
the  behavior  of  the  metabolite  is  complicated  due  to  its  strong  coupling  effects  at 
this  field  strength.  We  show  simulations  illustrating  the  behavior  of  the  resonance 
with  data  measured  from  phantom  and  in  vivo.  Echo  times  of  approximately  80  ms 
-  100  ms  gave  maximum  peak  amplitude  and  energy  of  the  metabolite  using  PRESS 
excitation.  We  also  show  simulations  using  2D  J-resolved  single  voxel  spectroscopy, 
which  provides  robust  detection.  Phantom  and  in  vivo  data  are  presented  to  illustrate 
the  spectral  pattern  of  the  J-resolved  Citrate  metabolite. 


Key  words:  citrate;  2D  J-resolved  spectroscopy,  prostate  cancer. 


1  Introduction 


Citrate  is  an  important  metabolite  often  used  in  Magnetic  Resonance  Spectroscopy 
(MRS)  and  MR  Spectroscopic  Imaging  (MRSI)  to  aid  the  detection  of  Prostate  Can¬ 
cer  (PCa)  along  with  the  Choline  and  Creatine  metabolites  (1).  In  addition  to  the 
morphological  information  provided  by  MR  Imaging  (MRI),  the  additional  informa¬ 
tion  gained  using  MRS/MRSI  increases  the  specificity  of  the  examination.  In  these 
exams,  the  level  of  (Choline  +  Creatine)/Citrate  is  regarded  as  a  marker  for  PCa 
(2).  To  date,  MRSI  protocols  for  PCa  detection  have  been  well  established  at  main 
magnetic  field  strength  of  1.5  Tesla  (T)  (3). 

The  advent  of  higher  field  strength  scanners  provide  many  new  extensions  from 
previous  1.5T  systems  due  to  the  inherent  increase  in  the  signal  to  noise  ratio  (SNR). 
For  PCa  exams  using  MRS/MRSI  methods,  this  advantage  can  be  exploited  in  var¬ 
ious  forms,  which  include  using  higher  spatial  resolution  acquisitions  to  increase  the 
accuracy  of  localization  of  the  cancerous  tissues.  Scan  times  can  also  be  made  shorter 
compared  to  1.5T  for  the  same  SNR  thereby  reducing  the  overall  MR  exam  time. 
The  extension  of  1.5T  MRSI  protocols  for  usage  in  3T  PCa  can  therefore  provide 
potential  merits. 

However,  the  process  of  advancing  to  higher  field  strength  accompanies  additional 
technical  considerations.  For  clinical  prostate  examinations  using  spectroscopic  tech¬ 
niques,  one  of  the  major  issues  that  arise  involves  the  detection  of  the  Citrate  reso¬ 
nance.  While  coupling  effects  can  be  largely  neglected  at  1.5T,  the  dominance  of  the 
strong  coupling  for  the  Citrate  resonance  at  3T  needs  consideration  in  selecting  the 
acquisition  parameters  if  reliable  detection  of  the  metabolite  is  desired.  At  3T,  the 
J-coupling  constant  is  closely  tied  with  the  chemical  shift  constant  of  the  metabo¬ 
lite  thereby  causing  dependance  of  the  spectra  as  a  function  of  echo  time  (TE)  (4). 
Constraints  in  the  radio  frequency  (RF)  field  can  also  exist  when  using  body  coil 
excitation  thereby  further  compromising  the  data  acquisition.  Therefore,  close  inves¬ 
tigation  of  the  behavior  of  the  Citrate  resonance  needs  to  be  studied  as  well  as  careful 
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considerations  in  selecting  the  image  protocols. 

This  study  involves  the  study  of  different  spectral  patterns  associated  with  the 
Citrate  resonance  at  3T.  We  show  the  dependance  as  a  function  of  the  echo  time 
using  simulations  based  on  the  full  density  matrix  calculation  and  using  phantom 
measurements  for  comparison.  Single  voxel  in  vivo  spectroscopic  data  are  acquired 
using  echo  times  selected  from  simulations  that  maximizes  the  detection  criteria.  Fur¬ 
thermore,  we  explore  the  use  of  2D  J-resolved  single  voxel  spectroscopic  sequence  (5) 
to  provide  robust  detection  of  the  Citrate  metabolite  and  to  illustrate  the  characteris¬ 
tics  of  the  J-resolved  spectral  pattern.  Simulations  and  measurements  are  performed 
with  phantom  while  in  vivo  data  are  collected  using  the  J-resolved  sequence.  This 
study  is  intended  to  provide  a  guideline  in  selecting  pulse  sequence  protocols  for  MR 
spectroscopy  to  aid  the  diagnosis  of  PCa  detection  at  3T. 


2  Methods 

Simulations  were  performed  for  the  Citrate  metabolite  by  solving  the  full  density 
matrix  of  strongly  coupled  two-spin  system  (6).  Different  echo  time  acquisitions 
were  simulated  assuming  a  PRESS  excitation  scheme.  The  timing  of  the  acquisition 
of  the  PRESS  sequence  was  assumed  to  be  90o-[fi2]-180o-[TE/2]-180o-[TE/2-f  125- 
acquisition,  where  TE  was  incremented  to  encode  echo  timing  dependance  (7).  For  our 
PRESS  sequence,  we  prescribed  tn  to  be  10  ms.  We  assumed  the  J-coupling  constant 
to  be  15.4  Hz  with  chemical  shift  of  0.12  ppm  (=16.6  Hz  at  3T)  for  simulations  (8). 
The  T2  value  was  assumed  to  be  200  ms  with  a  line  width  of  10  Hz.  The  simulations 
were  conducted  for  echo  times  starting  from  30  ms  and  ending  at  400  ms  with  10  ms 
intervals. 

Quantification  of  the  simulations  were  conducted  using  two  different  quantification 
metrics.  Firstly,  the  peak  amplitude  at  the  center  of  the  spectrum  (2.6  ppm  peak)  was 
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calculated  as  a  function  of  echo  time.  Secondly,  the  energy  was  calculated  by  inte¬ 
grating  the  square  of  the  spectrum.  Another  set  of  simulations  were  performed  which 
relates  to  the  power  constraints  in  the  RF  excitation  encountered  in  in  vivo  cases.  In 
such  situations,  reduced  flip  angle  techniques  are  often  required  (9).  Therefore,  we 
quantified  the  simulations  for  the  case  of  using  90°-180°-180°  flip  angle  excitation  as 
well  as  for  the  case  of  using  a  90°-167°-167°  excitation  often  used  in  in  vivo  exami¬ 
nations.  Experimental  data  were  acquired  from  a  Citrate  phantom  using  the  PRESS 
sequence  with  the  same  time  settings.  All  phantom  and  in  vivo  measurements  were 
also  made  with  this  reduced  flip  angle.  Data  from  different  echo  times  were  acquired 
with  a  3  second  repetition  time  (TR)  and  were  similarly  quantified  to  compare  with 
the  simulations.  All  data  were  acquired  on  a  3T  Signa  scanner  (GE  Medical  Systems, 
Waukesha,  WI). 

In  vivo  data  from  a  subject  suspicious  of  PCa  were  acquired  using  the  echo  time 
that  indicated  the  maximum  amplitude  and  energy  of  the  Citrate  metabolite  from 
the  simulation  and  phantom  measurement  results.  From  the  quantification  results, 
this  echo  time  for  maximum  detection  was  in  the  range  of  80-100  ms.  Echo  times 
below  50  ms  were  discarded  due  to  reasons  given  in  the  discussion  section.  Therefore 
in  vivo  single  voxel  spectroscopy  data  were  acquired  with  an  echo  time  of  90  ms.  Data 
were  acquired  for  8  minutes  from  a  voxel  of  size  1.5  cc  for  adequate  SNR.  All  in  vivo 
data  were  acquired  using  the  body  coil  for  excitation  followed  by  an  endorectal  coil 
signal  reception  (10).  All  in  vivo  studies  were  conducted  under  IRB  guidelines  and 
with  informed  consent. 

Simulations  and  phantom  measurements  were  also  conducted  for  a  2D  J-resolved 
acquisition.  In  both  cases,  the  echo  time  spacing  were  adjusted  to  be  7.8  ms  for  a 
total  of  64  steps  from  35  ms  to  534  ms  for  finer  resolution  in  the  FI  domain.  This 
resulted  in  a  2  Hz  spectral  resolution  with  a  bandwidth  of  128  Hz  in  the  FI  domain 
The  spectral  bandwidth  for  F2  was  5000  Hz  with  2048  data  point  acquisitions.  For 
the  actual  2D  J-resolved  measurements,  a  phantom  composed  of  Citrate,  Creatine, 
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and  Choline  metabolites  were  used  to  emulate  the  existence  of  cancerous  tissue. 

Finally,  in  vivo  data  were  collected  from  a  subject  suspicious  of  PCa  using  the 
2D  J-resolved  technique.  Single  voxel  2D  J-resolved  spectroscopic  data  were  acquired 
from  two  different  regions  at  the  peripheral  zone  of  the  prostate.  The  parameters 
were  maintained  as  with  the  simulations.  The  voxel  size  chosen  was  1cm  x  1.12cm 
x  1.08  cm  =  1.2  cc.  Four  acquisitions  were  averaged  per  echo  time  for  a  total  scan 
time  of  8  minutes  (TR  =  2  sec)  for  each  voxel. 


3  Results 

Figure  1  shows  the  Citrate  resonance  behavior  at  different  echo  times  obtained  from 
the  PRESS  excitation  simulation.  Unlike  1.5T,  due  to  the  strong  coupling  effects, 
very  different  resonant  behavior  can  be  observed  at  3T  as  a  function  of  echo  time. 
The  center  of  the  spectrum  can  be  seen  to  have  a  peak  at  very  short  echo  time  of  30 
ms,  a  negative  peak  near  90  ms  echo  time,  and  returning  to  a  positive  peak  near  270 
ms  but  at  a  reduced  amplitude  due  to  the  T2. 

Figure  2  quantifies  the  simulation  and  measurement  results  by  calculating  the 
relative  peak  value  at  the  2.6  ppm  location  (left)  and  the  energy  of  the  metabolite 
(right).  For  each  quantification  metric,  simulations  assuming  180°  flip  angle  (circled 
line),  simulations  assuming  reduced  flip  angle  excitation  (squared  line),  and  phantom 
measurement  results  (crossed  line)  are  shown.  Line  characteristics  of  the  measured 
data  obtained  from  the  phantom  are  seen  to  closely  resemble  the  results  obtained 
from  the  simulations  although  the  dynamic  range  is  smaller,  which  could  be  due  to 
inhomogeneities  or  other  inaccurate  modelling  of  the  metabolite.  In  both  quantifica¬ 
tions  however,  a  local  peak  can  be  obtained  with  echo  time  of  approximately  80  ms 
-  100  ms. 

In  Fig.  3,  in  vivo  single  voxel  measurement  data  acquired  with  90  ms  echo  time  is 
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presented.  The  excited  region  is  shown  on  the  T2  weighted  anatomical  image  (left) 
with  the  reconstructed  spectrum.  The  presence  of  Creatine  and  Choline  resonances 
can  be  seen  along  with  the  clearly  visible  inverted  peak  of  the  Citrate  metabolite  near 
the  2.6  ppm  region. 

Figure  4  shows  the  simulated  2D  J-resolved  Citrate  spectra  (right)  along  with  2D 
J-resolved  data  acquired  with  a  phantom  (left)  containing  Citrate  as  well  as  Choline 
and  Creatine  metabolites  for  different  FI  lines.  In  both  cases,  due  to  the  modulations 
occurring  as  a  function  of  echo  time,  resonances  are  clearly  seen  beyond  the  J(0)  line. 
The  spectra  for  both  cases  have  very  similar  spectral  patterns.  In  both  cases,  the 
J(0)  line  has  a  slight  negative  peak  at  the  Citrate  position  which  is  due  to  the  strong 
negative  peaks  at  echo  times  of  60  -  120  ms. 

Finally,  two  single  voxel  2D  J-resolved  spectra  from  in  vivo  subject  are  given  in 
Fig.  5.  The  two  regions  that  were  selected  are  shown  in  the  anatomical  T2  weighted 
image  given  along  with  the  2D  J-resolved  spectra.  The  spectra  obtained  from  the 
right  side  of  the  subject  (a)  displays  negligible  Citrate  metabolite  intensity  compared 
with  Creatine  and  Choline  resonances.  Compared  to  this  spectra,  the  2D  J-resolved 
spectra  from  the  left  side  of  the  subject  (b)  reveals  the  presence  of  Citrate  along 
with  the  Creatine  and  Choline  metabolites.  The  Citrate  metabolite  can  be  seen 
more  visibly  from  the  modulations  occurring  in  the  FI  lines.  These  two  comparisons 
show  that  with  the  2D  J-resolved  acquisition  method,  the  strongly  coupled  Citrate 
metabolite  can  be  resolved  while  the  presence  of  the  metabolite  can  be  confirmed 
robustly. 


4  Discussion 

The  successful  application  of  PCa  detection  using  MRS/MRSI  at  3T  depends  on 
the  proper  choice  of  pulse  sequence  parameters  based  on  the  knowledge  of  the  Citrate 
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resonance  behavior.  As  seen  from  Figs.  1  and  2,  to  maximize  the  resonance  amplitude 
and/or  energy,  good  choices  of  TE  would  be  below  35  ms  or  in  the  range  of  80  -  100 
ms.  In  the  former  case,  practical  issues  pertaining  to  minimizing  lipid  contamination 
using  saturation  pulses  and  eliminating  macromolecular  contributions  need  to  be 
established.  Also,  hardware  instabilities  at  short  echo  times,  for  example  gradient 
modulations  or  eddy  current  errors  can  deteriorate  the  quality  of  the  data  acquired. 
Meanwhile,  at  echo  times  between  80  -  100  ms,  lipid  contamination  can  be  reduced 
while  hardware  instabilities  can  be  less  problematic.  Therefore,  in  the  case  of  using 
a  fixed  echo  time  for  spectroscopy,  this  range  of  echo  times  can  be  advantageous. 

We  used  two  different  metrics  in  quantifying  the  Citrate  resonance.  Namely,  the 
peak  amplitude  at  the  center  of  the  spectrum  and  the  energy  of  the  metabolite, 
defined  by  the  integral  of  its  square  over  the  spectral  width.  Other  different  methods 
of  criteria  for  evaluating  the  outcome  can  also  be  established.  For  example,  in  the  case 
of  poor  B0  homogeneity,  line  broadening  can  cause  spectral  overlap.  Calculating  the 
energy  can  be  less  than  helpful  in  this  situation  while  negative  peaks  can  be  obscured. 
Rather,  obtaining  the  energy  of  the  main  central  lobe  can  give  more  robust  results  in 
such  cases. 

The  same  criteria  used  here  can  result  in  different  optimal  echo  times  depending  on 
the  pulse  sequence  parameters.  For  example,  simulations  assuming  a  90°-[TE/4]-180°- 
[TE/2]-180°-[TE/4]-acq  PRESS  sequence  resulted  in  maximum  (negative)  amplitude 
and  energy  at  a  slightly  shorter  echo  time  near  the  75  ms  range.  More  broadly 
speaking,  different  spectroscopic  sequences  such  as  J-PRESS  (4)  or  STEAM  (11)  can 
require  different  parameter  setting.  Finally,  restrictions  in  the  RF  field  can  also  yield 
different  results  which  is  a  nontrivial  practical  issue.  Even  with  the  simulations  at 
reduced  flip  angles  to  address  peak  limitation  issues,  the  B\  inhomogeneity  can  still 
be  in  effect  causing  parameter  settings  to  be  suboptimal.  Methods  to  reduce  RF 
inhomogeneity  and  power  constraints  is  highly  desirable  (12,  13). 

An  alternate  and  more  robust  method  for  detection  of  Citrate  is  to  use  J-resolved 
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spectroscopic  sequences  as  we  have  demonstrated  here.  In  this  case,  the  presence  of 
the  citrate  metabolite  is  well  demonstrated  at  lines  beyond  J(0),  which  provides  ro¬ 
bustness  in  the  estimation  of  the  presence  of  the  metabolite.  Other  metabolites  used 
as  markers  of  PCa  can  be  observed  in  the  J(0)  line.  In  addition  to  the  benefit  of  resolv¬ 
ing  the  J-coupling  for  the  detection  of  Citrate,  it  has  been  shown  that  J-resolved  spec¬ 
troscopy  can  further  differentiate  other  metabolites,  for  example,  polyamines  which 
co-resonate  near  the  Creatine  and  Choline  peaks  (14,  15).  This  eventually  leads 
to  additional  physiologic  information  concerning  the  prostate  tissue.  Furthermore, 
J-resolved  spectroscopy  has  been  used  to  reduce  sidebands  artifacts  (16,  17).  The  ac¬ 
quisitions  from  multiple  echo  times  can  also  help  determine  the  T2  of  the  metabolites 
of  interest  as  well  as  the  water  content(14). 

For  application  of  J-resolved  techniques  in  the  clinical  settings,  a  means  of  quan¬ 
tifying  the  amount  of  the  Citrate  metabolite  needs  to  be  established.  Due  to  the 
modulation  of  the  citrate  resonance  shape  as  a  function  of  echo  time,  a  straight  for¬ 
ward  method  of  quantifying  the  amount  is  not  available.  Rather,  the  relative  amount 
can  be  inferred  from  the  peaks  at  J(±4),  J(±8),  or  other  lines  in  FI  using  a  model 
based  approach. 

A  clinically  more  useful  examination  using  J-resolved  spectroscopy  would  be  to 
have  information  regarding  the  spatial  distribution  of  the  metabolites.  Using  straight 
forward  phase  encoding  methods  for  gathering  spatial  information  can  lead  to  in¬ 
creased  scan  time  not  suitable  for  in  vivo  studies.  To  overcome  these  issues,  clever 
choice  of  pulse  sequence  parameters  needs  to  be  made.  For  example,  a  multiple  spin 
echo  sequence  can  be  used  to  collect  additional  echo  time  data  from  one  excitation 
pulse.  Also,  since  the  Citrate  resonance  is  mostly  concentrated  near  the  center  of 
the  J(0)  line,  adjusting  the  sampling  can  be  used  to  effectively  decrease  the  scan 
without  deteriorating  the  spectra  by  using  undersampling  or  variable  rate  sampling 
techniques  (18).  Another  alternative  method  to  increasing  spatial  information  is  to 
use  time  varying  readout  gradients  (19,  20).  It  has  been  shown  that  using  read- 
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out  gradients  can  dramatically  decrease  the  minimum  total  scan  time.  This  can  be 
effectively  used  to  gain  the  additional  information  needed  for  J-resolved  techniques. 

5  Conclusion 

An  analysis  of  the  Citrate  metabolite,  which  is  useful  for  the  detection  of  Prostate 
cancer,  using  full  density  matrix  has  been  conducted  for  applications  at  3T.  Phantom 
results  axe  provided  for  comparison.  Using  peak  estimate  and  energy  calculations, 
echo  time  in  the  range  of  80  -  100  ms  gave  maximum  results  suggesting  single  voxel 
or  chemical  shift  imaging  experiments  using  PRESS  excitation  to  be  conducted  at 
this  echo  time.  Also,  we  have  shown  that  2D  J-resolved  spectroscopy  provides  an 
elegant  method  to  detect  the  Citrate  metabolite  effectively.  Strong  coupling  which 
leads  to  modulations  as  a  function  of  echo  time  is  easily  observed  using  2D  J-resolved 
spectroscopic  sequences.  In  vivo  measurements  confirms  the  observation  of  Citrate 
with  patterns  similar  to  simulation  and  phantom  studies. 
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Table  and  Figure  Captions 


Figure  1.  Citrate  metabolite  characteristics  for  different  echo  times  at  3T.  The  Cit¬ 
rate  resonance  behavior  was  simulated  assuming  a  PRESS  excitation  with  J-coupling 
constant  of  15.4  Hz  and  chemical  shift  value  of  0.12  ppm  (16.6  Hz  at  3T).  The  behav¬ 
ior  is  very  different  from  1.5  T  due  to  the  strong  coupling  effects.  The  central  lobe 
modulates  as  a  function  of  echo  time,  with  positive  peak  at  short  TE«30  ms  followed 
by  a  negative  peak  near  the  90  ms  range.  Another  positive  peak  can  be  observed  at 
270  ms  range  but  decreased  in  intensity  due  to  T2  effects,  which  was  assumed  to  be 
200  ms. 

Figure  2.  Quantitative  results  obtained  from  simulations  and  measurement  from 
phantom  as  a  function  of  echo  time.  Graphs  are  plotted  for  simulated  results  assuming 
90°-180o-180°  excitation  (o),  simulated  results  assuming  reduced  flip  angle  90°-167°- 
167°  excitation  (o),  and  measured  data  acquired  from  phantom  (x).  Left)  value  of 
center  peak  of  the  spectrum  at  2.6  ppm.  Right)  energy  of  the  Citrate  spectrum. 
The  comparison  in  both  cases  shows  good  agreement  between  the  simulations  and 
phantom  measurements.  Also,  it  provides  a  good  guideline  for  selecting  the  echo 
time  for  Citrate  metabolite  detection.  In  both  cases,  echo  times  near  80  -  100  ms 
range  have  a  local  peak. 

Figure  3.  In  vivo  single  voxel  spectroscopy  measurement.  Measurement  region  was 
selected  from  the  T2  weighted  anatomic  image  as  shown  using  the  PRESS  sequence 
at  an  echo  time  of  90  ms  where  an  inverted  peak  of  the  Citrate  metabolite  should 
be  observed.  As  the  result  spectrum  suggests,  the  presence  of  the  Citrate  metabolite 
can  be  seen  at  the  2.6  ppm  region  as  an  inverted  peak.  Other  metabolites  of  interest 
for  detection  of  PCa,  namely  the  Creatine  and  Choline  peaks  can  also  be  observed. 

Figure  4.  2D  J-resolved  spectra  obtained  from  simulations  (right)  and  phantom  mea¬ 
surements  (left)  at  3T  using  the  PRESS  sequence.  The  echo  time  interval  was  7.8 
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ms  starting  from  35  ms  for  64  steps.  For  the  phantom  measurement,  Choline  and 
Creatine  metabolites  were  added  to  mimic  PCa.  Due  to  modulations  as  a  function 
of  echo  time,  the  J-resolved  spectra  show  strong  signal  of  the  Citrate  metabolite  at 
lines  outside  of  J(0).  Simulated  J-resolved  spectra  and  phantom  measured  J-resolved 
spectra  are  shown  with  similar  spectral  patterns.  In  this  respect,  the  detection  of  the 
Citrate  resonance  can  be  made  robustly  using  J-resolved  acquisitions. 

Figure  5.  Single  voxel  2D  J-resolved  spectroscopy  results  obtained  from  an  in  vivo 
subject  suspicious  of  PCa.  Two  voxels  were  selected  for  the  single  voxel  examinations 
as  shown  in  the  T2  weighted  anatomical  images.  The  corresponding  J-resolved  spectra 
are  shown  on  the  bottom  of  each  image.  In  (a),  even  though  the  presence  of  Creatine 
and  Choline  metabolites  are  evident,  there  is  no  visible  Citrate.  As  for  the  region 
shown  in  (b),  the  Citrate  resonance  is  visible  (2.6  ppm  region  from  J(-10)  to  J(10) 
line)  in  the  J-resolved  spectra  while  other  metabolites  are  also  present.  This  shows 
the  J-resolved  spectroscopy  can  be  useful  for  the  detection  of  the  Citrate  metabolite 
as  well  as  other  metabolites  of  interest  in  the  detection  of  PCa. 
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Figure  1:  Citrate  metabolite  characteristics  for  different  echo  times  at  3T.  The  Cit¬ 
rate  resonance  behavior  was  simulated  assuming  a  PRESS  excitation  with  J-coupling 
constant  of  15.4  Hz  and  chemical  shift  value  of  0.12  ppm  (16.6  Hz  at  3T).  The  behav¬ 
ior  is  very  different  from  1.5  T  due  to  the  strong  coupling  effects.  The  central  lobe 
modulates  as  a  function  of  echo  time,  with  positive  peak  at  short  TE«30  ms  followed 
by  a  negative  peak  near  the  90  ms  range.  Another  positive  peak  can  be  observed  at 
270  ms  range  but  decreased  in  intensity  due  to  T2  effects,  which  was  assumed  to  be 
200  ms. 
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Figure  2:  Quantitative  results  obtained  from  simulations  and  measurement  from 
phantom  as  a  function  of  echo  time.  Graphs  are  plotted  for  simulated  results  as¬ 
suming  90°-180°-180°  excitation  (o),  simulated  results  assuming  reduced  flip  angle 
90°-167°-167°  excitation  (o),  and  measured  data  acquired  from  phantom  (x).  Left) 
value  of  center  peak  of  the  spectrum  at  2.6  ppm.  Right)  energy  of  the  Citrate  spec¬ 
trum.  The  comparison  in  both  cases  shows  good  agreement  between  the  simulations 
and  phantom  measurements.  Also,  it  provides  a  good  guideline  for  selecting  the  echo 
time  for  Citrate  metabolite  detection.  In  both  cases,  echo  times  near  80  —  100  ms 
range  have  a  local  peak. 
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Figure  3:  In  vivo  single  voxel  spectroscopy  measurement.  Measurement  region  was 
selected  from  the  T2  weighted  anatomic  image  as  shown  using  the  PRESS  sequence 
at  an  echo  time  of  90  ms  where  an  inverted  peak  of  the  Citrate  metabolite  should 
be  observed.  As  the  result  spectrum  suggests,  the  presence  of  the  Citrate  metabolite 
can  be  seen  at  the  2.6  ppm  region  as  an  inverted  peak.  Other  metabolites  of  interest 
for  detection  of  PCa,  namely  the  Creatine  and  Choline  peaks  can  also  be  observed. 


Figure  4:  2D  J-resolved  spectra  obtained  from  simulations  (right)  and  phantom  mea¬ 
surements  (left)  at  3T  using  the  PRESS  sequence.  The  echo  time  interval  was  7.8 
ms  starting  from  35  ms  for  64  steps.  For  the  phantom  measurement,  Choline  and 
Creatine  metabolites  were  added  to  mimic  PCa.  Due  to  modulations  as  a  function 
of  echo  time,  the  J-resolved  spectra  show  strong  signal  of  the  Citrate  metabolite  at 
lines  outside  of  J(0).  Simulated  J-resolved  spectra  and  phantom  measured  J-resolved 
spectra  are  shown  with  similar  spectral  patterns.  In  this  respect,  the  detection  of  the 
Citrate  resonance  can  be  made  robustly  using  J-resolved  acquisitions. 
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Figure  5:  Single  voxel  2D  J-resolved  spectroscopy  results  obtained  from  an  in  vivo 
subject  suspicious  of  PCa.  Two  voxels  were  selected  for  the  single  voxel  examinations 
as  shown  in  the  T2  weighted  anatomical  images.  The  corresponding  J-resolved  spectra 
are  shown  on  the  bottom  of  each  image.  In  (a),  even  though  the  presence  of  Creatine 
and  Choline  metabolites  are  evident,  there  is  no  visible  Citrate.  As  for  the  region 
shown  in  (b),  the  Citrate  resonance  is  visible  (2.6  ppm  region  from  J(-10)  to  J(10) 
line)  in  the  J-resolved  spectra  while  other  metabolites  are  also  present.  This  shows 
the  J-resolved  spectroscopy  can  be  useful  for  the  detection  of  the  Citrate  metabolite 
as  well  as  other  metabolites  of  interest  in  the  detection  of  PCa. 
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1  Introduction 


1.1  Molecular  and  functional  imaging 

For  much  of  the  last  century,  medical  imaging  has  been  focused  on  faster  and  more 
detailed  anatomic  pictures  of  the  human  body.  The  accomplishment  of  the  visible 
human  project  of  the  National  Library  of  Medicine  (http://www.nlm.nih. gov/ 
research/visiblel  represents  perhaps  one  of  the  most  important  milestones  in  these 
developments.  With  the  goal  of  producing  a  system  of  knowledge  structures  that 
transparently  links  visual  knowledge  forms  to  symbolic  knowledge  formats  such  as 
the  names  of  body  parts,  a  complete,  anatomically  detailed,  3D  representations  of 
the  normal  male  and  female  human  bodies  were  rendered  based  on  transverse  CT, 
MR  and  cryosection  images  of  male  and  female  cadavers.  Medical  imaging  has 
been  an  integral  part  of  radiation  therapy  since  the  discovery  of  x-rays  and  the 
imaging  techniques,  such  as  X-ray,  CT,  MR1  and  ultrasound  imaging,  are  the 
foundation  for  the  modem  radiation  therapy  modalities  that  are  routinely  used  in  the 
clinics,  such  as  3D  conformal  radiation  therapy,  intensity  modulated  radiation 
therapy  (IMRT),  stereotactic  radiosurgery,  and  brachytherapy.  Indeed,  the 
development  of  radiation  therapy  has  strongly  relied  on  the  imaging  technolgy  and, 
historically,  almost  every  major  advancement  in  imaging  science  would  bring 
radiation  therapy  to  a  new  level. 

In  general,  medical  imaging  is  involved  in  all  key  steps  of  radiation  treatment 
(figure  1).  One  of  the  most  important  uses  of  imaging  techniques  is  the  delineation 
of  a  tumor  target.  Despite  the  tremendous  successes,  the  anatomic  imaging 
techniques  such  as  CT/MRI/US  are  inherently  deficient  in  that  they  can  only  reveal 
spatial  changes  in  physical  properties  and  fail  to  provide  basic  biological 
information  that  is  much  needed  for  the  optimal  management  of  the  patients. 
Clinically,  tumor  biology  plays  an  important  role  in  the  diagnosis,  treatment 
decision-making,  and  assessment  of  therapeutic  response  of  various  diseases.  It  is 
thus  highly  desiable  to  develop  imaging  techniques  capable  of  revealing  the  spatial 
biology  distribution  of  the  patients.  Toward  this  goal,  a  new  brach  of  science, 
referred  to  as  molecular  imaging,  is  emergying  as  a  result  of  research  efforts  in 
cellular  biology  and  imaging  techniques  over  the  years.  The  development  of 
cellualar  and  molecular  imaging  provides  significant  opportunities  for  radiation 
discipline  to  take  patient’s  biological  information  into  radiation  therapy  treatment 
decision-making  process  and  to  truly  individualize  the  cancer  radiotherapy. 
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1 .2.  IMRT  as  a  means  of  producing  biologically  conformal  dose  distributions 
IMRT  is  an  advanced  form  of  external  beam  irradiation  and  represents  a  radical 

change  in  radiation  oncology  practice  This  new  process  of  treatment  planning 
and  delivery  shows  significant  potential  for  improving  therapeutic  ratio  and  offers  a 
valuable  tool  for  dose  escalation  and/or  radiation  toxicity  reduction.  Preliminary 
published  results  and  unpublished  results  from  several  institutions  indicate  that  with 
IMRT,  radiation  doses  to  sensitive  structures  can  be  reduced  significantly  while 
maintaining  adequate  target  dose  coverage  4-15  Because  in  many  clinical  situations 
the  dose  to  the  tumor  volume  is  limited  by  the  tolerance  doses  of  the  sensitive 
structures,  it  is  considered  likely  that  IMRT  will  improve  local  control  and  lead  to  an 
increase  in  survival  rate  for  certain  cases  through  dose  escalation.  In  addition,  IMRT 
has  the  potential  to  improve  the  efficacy  of  treatment  planning  and  delivery  in 
routine  clinical  practice  with  the  use  of  computerized  planning  and  treatment 
process.  For  details  about  IMRT  inverse  treatment  planning,  deivery  and  quality 
assurance,  we  refer  the  readers  to  the  related  chapters  of  this  book. 

In  IMRT,  each  incident  beam  is  divided  into  a  number  of  beamlets  (typically, 
the  size  of  beamlet  is  lcmXlcm),  allowing  us  to  modify  the  dose  distribution  on  an 
individual  beamlet  level.  Using  IMRT,  it  is  possible  to  produce  not  only  spatially 
uniform  but  also  non-uniform  dose  distributions.  Recently,  Ling  et  al  and  several 
other  researchers  16-21  have  emphasized  the  technical  capability  of  "dose  painting" 
and  "dose  sculpting"  offered  by  IMRT,  which  allows  customized  dose  delivery  to  the 
target  volume(s)  with  centimeter  or  even  sub-centimeter  spatial  resolution.  Using 
functional  and  molecular  imaging  techniques  to  identify  spatial  metabolic 
distribution  and  hence  guide  the  delivery  of  radiation  represents  a  paradigm  shift  in 
radiation  oncology  and  this  type  of  “biologically”  conformal  radiation  therapy  may 
provide  a  significant  opportunity  to  improve  conventional  IMRT  treatment.  A  timely 
question  is  how  to  integrate  the  state-of-the-art  functional  imaging  technologies  into 
radiation  therapy  techniques  such  as  IMRT  to  positively  impact  clinical  cancer 
management.  The  purpose  of  this  Chapter  is  to  review  recent  progress  in  this 
endeavor  and  identify  the  important  issues  in  the  development  of  biologically 
conformal  radiation  therapy. 
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2.  Biologically  conformal  IMRT 

Current  IMRT  treatment  plan  optimization  is  based  on  the  assumption  of  uniform 
biology  distribution  within  the  target  volume  and  is  aimed  at  achieving 
geometrically  conformal  dose  distributions  under  the  guidance  of  CT/MRI  images. 
In  reality,  it  has  long  been  recognized  that  the  spatial  distribution  of  biological 
properties  in  most  tumors  and  normal  tissues  are  heterogeneous.  With  the  advent  of 
various  molecular  and  functional  imaging  techniques,  it  is  now  possible  to  map  out 
biology  distribution  on  a  patient  specific  basis.  To  use  the  spatially  heterogeneous 
biology  information  derived  from  the  new  imaging  modalities  to  guide  IMRT  dose 
painting  and  sculpting  process,  several  key  problems  need  to  be  resolved.  In  general, 
the  molecular/functional  imaging-guided  IMRT  generally  favors  non-uniform  dose 
distributions  and  requires  a  plan  optimization  formalism  in  voxel  domain  to  deal 
with  the  biological  heterogeneity.  In  addition,  new  method  of  specifying  the  desired 
doses  and  a  mechanism  for  inter-  and  intra-structural  tradeoff  must  be  introduced  to 
efficiently  produce  metabolically/functionally  conformal  doses.  In  figure  2  we  list  the 
general  steps  of  biologically  conformal  IMRT  treatment.  Each  of  the  steps  in  figure  2 
is  discussed. 

2.1  Integration  of  functional  and  molecular  imaging  into  IMRT  planning 

The  area  of  molecular/functional  imaging  is  rapidly  evolving  22-24.  Many  of  the 
molecular  imaging  modalities  (such  as  fluorescent  and  bioluminescent  imaging, 
optical  imaging,  SPECT/PET  with  novel  isotopes/contrast  agents  targeting  some 
specific  molecular  markers,  MR  spectroscopic  imaging  (MRSI))  are  being 
developed  for  tumour  specific  imaging  and  deployed  into  clinical  practice. 
Presently,  MRSI,  PET/SPECT  and  micro-bubble  based  ultrasound  are  perhaps  the 
most  mature  modalities  and  available  for  guiding  radiation  therapy  treatment. 
Details  on  various  molecular  and  functional  imaging  modalities  have  been  given  in 
Chapter  3  and  will  not  be  repeated  here.  We  refer  the  readers  to  that  Chapter  and 
references  therein.  The  remaining  of  this  Chapter  will  be  focused  on  the  issues 
related  to  the  integration  of  the  new  imaging  modalities  into  radiation  treatment 
planning. 

2.2  Image  registration 

Radiation  therapy  treatment  planning  is  mainly  CT  image-based  because  it 

provides  complete  geometric  data  and  electron  density  information  for  accurate  dose 
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calculation.  To  utilize  the  biological  information  derived  from  the  new  image 
modalities,  we  must  map  the  imaging  data  onto  treatment  planning  CT  images.  The 
level  of  complexity  of  image  coregistration  depends  on  the  imaging  techniques 
involved  and  specific  software  tools  often  need  to  be  developed  in  order  to  use  some 
of  the  new  imaging  modalities,  such  as  fluorescent  images,  endoscopic  images  and 
endorectal  images.  Sometimes,  deformable  model-based  image  registration  is 
required  if  the  shape(s)  of  the  involved  organs  are  deformed  from  its  normal  shape. 

Let  us  take  endorectal  MRSI  as  an  example.  The  introduction  of  endorectal 
surface  coils  significantly  improve  spatial  resolution  and  signal-to  noise  ratio  (SNR) 
of  prostate  MR  imaging  and  allows  evaluation  of  the  tumor  location,  tumor  volume, 
capsular  penetration,  invasion  of  neurovascular  bundle,  and  seminal  vesicle 
involvement,  which  is  crucial  for  accurate  treatment  planning.  Endorectal-coil  based 
MRSI  has  also  been  shown  effective  in  distinguishing  between  areas  of  cancer  and 
normal  prostatic  epithelium  through  differences  in  [choline+  creatine]/citrate  ratio 

25-28  However,  the  use  of  endorectal  probe  inevitably  distorts  the  prostate  and 
other  soft  tissue  organs,  making  it  impossible  to  directly  fuse  the  acquired  image 
data  onto  treatment  planning  CT.  In  figure  3  we  show  the  difference  between 
endorectal  coil-based  MRI  defined  and  CT-defined  prostate  volume  29.  jn  order  to 
fuse  MRI/MRSI  with  treatment  planning  CT,  it  is  needed  to  develop  an  effective 
deformable  image  registration  procedure.  Otherwise,  the  gain  from  the  use  of  the 
state-or-the-art  imaging  techniques  may  be  lost  due  to  the  inferior  performance  of 
image  registration. 

Zaider  et  al  30  have  reported  a  translation  and  scaling  based  registration 
method  to  map  MRS  positive  volumes  onto  the  CT  and  ultrasound  images.  In  their 
approach,  the  coordinates  of  the  boundary  and  the  center  of  mass  were  used  to 
linearly  interpolate  the  positions  of  the  mapped  voxels.  A  larger  discrepancy  was 
found  for  regions  with  more  severe  distortion  (>  4mm).  Lian  et  al  29,  31  have 
developed  an  effective  deformable  image  registration  algorithm  to  map  the 
MRI/MRSI  information  obtained  using  a  rigid  or  inflatable  endorectal  probe  onto 
CT  images  and  to  verify  the  accuracy  of  the  registration  by  phantom  and  patient 
studies.  For  this  purpose,  a  thin  plate  spline  (TPS)  transformation  first  introduced  by 

Bookstein  32  was  implemented  to  establish  voxel-to-voxel  correspondence  between 
a  reference  image  and  a  floating  image  with  deformation.  The  idea  is  to  find  a 

continuous  transformation  within  a  suitable  Hilbert  space  to  minimize  the  landmark 
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difference  in  two  images.  The  detailed  description  of  the  TPS  transformation  can  be 
found  in  Bookstein’s  original  paper  32.  To  access  the  quality  of  the  registration,  an 
elastic  phantom  with  a  number  of  implanted  fiducial  markers  was  designed. 
Radiographic  images  of  the  phantom  were  obtained  before  and  after  a  series  of 
intentionally  introduced  distortions.  After  mapping  the  distorted  phantom  to  the 
original  one,  the  displacements  of  the  implanted  markers  were  measured  with 
respect  to  their  ideal  positions  and  the  mean  error  was  calculated.  Phantom  studies 
showed  that  using  the  deformable  registration  method  the  mean  landmark 
displacement  error  was  0.62  ±  0.39  mm  when  the  distortion  was  of  the  order  of 
23.07  mm.  A  deformable  model  seems  to  b  necessary  to  faithfully  map  the 
metabolic  information  onto  the  treatment  CT  images.  When  a  non-deformable 
method  based  on  a  rigid-body  transformation  and  scaling  was  used  for  the  same 
distortion,  the  mean  displacement  of  the  fiducials  with  respect  to  their  actual 
positions  was  found  to  be  as  large  as  12.95  ±  0.57  mm.  In  patient  studies,  CT  images 
of  two  prostate  patients  were  acquired,  followed  by  3  Tesla  (3T)  MR  images  with  a 
rigid  endorectal  coil.  For  both  patient  studies,  significantly  improved  registration 
accuracy  was  achieved.  The  prostate  centroid  position  displacement  was  0.58  ±0.10 
mm  and  the  coincidence  index  was  92.6  ±  5.1  %  when  a  TPS  transformation  was 
used.  Different  from  the  non-deformable  approach,  the  TPS-based  registration 
accommodates  the  organ  distortion  and  enables  us  to  achieve  significantly  higher 
MR/MRSI  and  CT  image  registration  accuracy.  More  advanced  finite  element 
method  is  also  developed  to  attack  the  problem^. 

2.3  Quality  assurance  of  molecular  and  functional  imaging  modalities 

Any  new  imaging  modality  requires  validation  and  quality  assurance  to  ensure 
that  the  obtained  images  faithfully  reflect  the  reality.  In  anatomical  imaging, 
surrogate  phantoms  have  been  widely  used  for  assessing  the  geometric  and  physical 
(e.g.,  electron  density)  properties  of  the  images.  For  radiation  therapy  application, 
Mutic  et  al  have  reported  a  simple  design  of  PET  phantom  to  validate  the  image 
registration  of  PET  and  CT  images34.  Generally  speaking,  for  a  biological  imaging 
modality,  validation  of  geometric  accuracy  represents  only  one  facet  of  the  problem. 
The  accuracy  of  the  pixel  values  of  the  imaging  modality  also  needs  our  attention. 
While  the  specific  meaning  of  the  pixel  values  depends  on  the  modality,  let  us  take 
an  MRSI  phantom  (figure  4)  constructed  by  Hunjan  et  al  as  an  example  to  illustrate 
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the  basic  idea.  The  multi-modality,  multi-purpose  phantoms  is  suitable  for  quality 
assurance  testing  of  fusion  data  from  MRI,  MRSI  and  CT  images35.  The  phantom 
contains  fiducial  markers  that  are  simultaneously  MR,  MRS,  and  CT-visible.  To 
examine  the  accuracy  of  MRSI  for  brain  tumor,  the  phantom  was  filled  with  a  brain- 
mimicking  solution  with  an  insert  holding  8  vials  containing  calibrated  solutions  of 
precisely  varying  metabolite  concentrations  that  emulated  increasing  grade/density 
of  brain  tumor.  Metabolite  ratios  calculated  from  fully  relaxed  ID,  2D  and  3D  MRS 
data  for  each  vial  were  compared  to  calibration  ratios  acquired  in  vitro  using  a  9.4 
Tesla  MR  spectrometer.  Figure  5  shows  an  axial  scout  scan  of  the  MRS  metabolite 
ratio  quantitation  standard  showing  the  calibration  vials  1-8.  The  resulting  single 
voxel  MR  spectra  are  shown  inset  next  to  corresponding  vials  and  a  linear  fit 
between  the  Choline/NAA  ratio  of  the  calibration  solutions  obtained  at  9.4  T  versus 
the  calibration-solution-filled  vials  inside  the  phantom  obtained  at  1.5  T.  For 
detailed  information  on  the  design  of  the  phantom  and  measurements,  please  refer 
35,  36, 

2.4  Inverse  planning 

In  general,  molecular/functional  imaging  could  impact  the  current  radiation 
therapy  treatment  in  two  fundamental  aspects  16,  20,  37,  First,  it  offers  an  effective 
means  for  us  to  more  accurately  delineate  the  tumor  and  better  define  the  treatment 
volume.  Secondly,  it  provides  valuable  spatial  metabolic  information  in  the  tumor 
and  sensitive  structures.  While  it  is  straightforward  to  modify  the  radiation  portals  to 
accommodate  any  changes  in  treatment  volume,  new  methods  of  dose  optimization 
and  medical  decision-making  must  be  developed  to  take  full  advantage  of  the 
metabolic  information  and  IMRT.  We  have  recently  introduce  the  concept  of  4D 
inverse  planning  and  demonstrated  the  integration  of  MRSI  into  IMRT  planning. 
The  goal  of  this  type  of  treatment  scheme  is  to  achieve  biologically  conformal  doses, 
instead  of  the  geometrically  conformal  dose  distribution  sought  by  conventional 
radiation  therapy  planning.  We  showed  that,  under  the  guidance  of  MRSI  metabolic 
maps,  it  is  possible  to  prescribe  a  higher  dose  where  there  is  resistance  and/or  where 
there  are  dense  tumor  cell  populations.  Similarly,  the  technique  also  allows  for 
differential  sparing  of  sensitive  structures  accounting  for  functionally  important 
regions.  A  few  important  issues  related  to  new  inverse  planning  scheme  are  outlined 
in  the  following. 
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2.4.1  Objective  function 

The  objective  function  is  used  in  inverse  planning  to  quantitatively  rank  the 
candidate  plans  and  the  optimization  of  the  function  yields  the  optimal  input  beam 
parameters.  Ideally,  an  objective  function  should  rank  a  given  solution  in  a  way 
consistent  with  the  clinical  judgment.  Despite  of  intense  research  effort  in  modeling 
the  clinical  decision-making  strategies  38,  39  40  41,  42;  appropriate  form  of  the 
objective  function  for  clinical  IMRT  planning  remains  illusive  even  for  conventional 
IMRT  inverse  planning.  Our  recent  study  has  indicated  that  the  IMRT  plans  derived 
from  the  existing  inverse  planning  techniques  are,  at  best,  sub-optimal.  For 
convenience,  it  seems  appropriate  to  classify  the  currently  available  dose 
optimization  methods  into  four  categories:  (i)  dose-based43-48-  (H)  clinical 
knowledge-based49;  (Hi)  EUD-feased  50-53-  ancj  (,v)  TCP/NTCP-based38,  54,  55, 
The  underlying  difference  between  these  models  lies  in  that  what  endpoint  are  used 
to  quantify  the  treatment  or  what  fundamental  quantities  are  used  to  define  the 
optimality  of  a  plan.  In  reality,  each  type  of  inverse  planning  formalism  has  pros  and 
cons.  The  dose-based,  EUD-based  and  TCP/NTCP-based  optimizations  have  been 
discussed  extensively  in  the  previous  Chapters.  We  focus  our  discussion  on  the  2nd 
category  here. 

To  illustrate  the  need  for  a  biologically  more  sensible  yet  clinically  more 
relevant  approach,  let  us  take  parotid  glands  as  an  example.  It  is  well  known  that  the 
clinical  end  point  is  the  same  if  the  glands  are  irradiated  15Gy  to  67%,  or  30Gy  to 

30%,  or  45Gy  to  24%  of  the  total  volume56.  A  quadratic  objective  function  or  alike 
would  rank  the  three  scenarios  completely  differently.  The  use  of  dose-volume 
constraints  does  not  change  the  ranking  because  they  act  as  “boundary  conditions” 
only  limit  the  search  space.  A  minimum  requirement  for  a  sensible  objective 
function  is  that  it  should  be  consistent  with  the  existing  clinical  outcome  data.  While 
the  biologically  based  models  are  more  relevant  for  radiotherapy  plan  ranking,  the 
dose-response  function  of  various  structures  is  not  sufficiently  understood  and,  at 
this  point,  there  is  considerable  controversy  about  the  models  for  computing  dose- 
response  indices  and  their  use  in  optimization.  In  reality,  it  is  difficult  for  clinicians 
to  specify  the  optimization  criteria  in  terms  of  certain  dose  response  indices  (e.g., 
TCP,  NTCP  and  P+).  This  issue  becomes  compounded  when  two  or  more 
independently  optimized  plans  are  to  be  combined.  Towards  developing  a  clinically 
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practical  inverse  planning  formalism,  we  have  proposed  a  new  type  of  the  inverse 
planning  formulation  with  integration  of  clinical  end  point  data  and  showed  that  the 
clinical  knowledge-based  model  more  objectively  ranks  radiation  therapy  plans  and 
makes  it  possible  to  obtain  truly  optimal  IMRT  plans  with  much  reduced  efforts. 

The  objective  function  is  based  on  the  concept  of  equivalent  volume57'59, 
for  an  individual  voxel  /', 

(A ,  (1) 

where  F  is  the  volume  of  the  voxel  i,  A  and  A*/ are  the  dose  and  reference  dose  at 
the  voxel,  respectively,  n  is  a  parameter  that  represents  the  biological  characteristic 
of  the  organ.  For  a  sensitive  structure,  n  is  a  small  positive  number  (0<n<l)  and  the 
value  of  parameter  n  reflects  the  architecture  (serial  or  parallel)  of  the  sensitive 
structure.  For  a  target,  n  should  be  assigned  with  a  small  negative  value  (-l<n<0).  In 
order  to  model  the  volumetric  effect,  the  objective  function  should  take  the  form  of 
/,=/,({(  AVeff),}).  (2) 

In  figures  6  and  7  we  compare  IMRT  plans  obtained  using  the  clinical  knowledge- 
based  and  conventional  quadratic  objective  functions.  As  seen  from  the  figures,  dose 
distribution  was  significantly  improved  as  a  result  of  the  more  adequate  modeling. 

2.4.2  Relation  between  metabolic  abnormality  level  and  radiation  dose 

An  important  task  in  biologically  conformable  radiation  therapy  is  to  quantify  the 
tumor  burden  and  relate  the  metric  to  the  radiation  dose.  In  Sec.  2.4.6  we  will  derive 
such  a  relation  based  on  radiobiology  model  under  the  assumption  that  all  necessary 
model  parameters  are  known.  At  present,  the  radiobiology  parameters  are  spares  and  one 
should  perhaps  take  a  less  precise  yet  more  practical  approach.  Generally  speaking,  the 
relation  between  the  abnormality  level  and  radiation  dose  can  in  principle  be  determined 
experimentally  or  through  analysis  of  animal  and  hypothesis-driven  clinical  data,  which 
is  similar  to  the  establishment  of  the  empirical  radiation  dose  prescriptions  for  different 
disease  sites  in  our  current  clinical  practice.  A  linear  relation  between  the  dose  and 
metabolic  abnormality  levels  20 

D?(n)=D0'+KM(n),  (3) 

was  assumed  in  our  previous  study,  where  Df(n)  is  the  prescribed  target  dose  at  voxel  n, 
M(n)  is  the  abnormality  level  at  the  voxel,  a:  is  an  empirical  coefficient,  Do  can  be 
regarded  as  the  conventional  prescription  dose  when  functional  imaging  information  is 


9 


not  available  or  when  the  abnormality  level  is  minimum,  M(n)=0.  The  bottom  line  is 
that  no  subvolume  in  the  tumor  should  receive  a  dose  less  than  conventionally  prescribed 
dose  Do  unless  it  is  clinically  justifiable.  For  a  given  organ,  we  postulated  that  the 
tolerance  dose  is  related  to  the  functional  importance  by 

Dc(n)=  Do-ocK(n) ,  (4) 

where  Dc(n)  is  the  tolerance  dose  at  voxel  n,  K(n)  is  the  functional  importance  at  the 
voxel,  a  is  an  empirical  coefficient.  Do  represents  the  tolerance  dose  corresponding  to 
the  situation  when  functional  distribution  information  is  not  available  or  when  the 
functional  importance  is  minimum,  K(n)=0. 

We  emphasize  that  the  above  two  relations  are  somewhat  ad  hoc  and  may  need  to 
be  refined  as  more  knowledge  are  gained.  For  treatment  of  dose-response  neoplasm, 
such  as  prostate  cancer,  it  seems  to  be  a  good  strategy  to  attempt  to  escalate  the  dose  to 
those  high  tumor  burden  points  as  high  as  possible  while  keeping  the  normal  tissue 
complications  below  a  certain  level.  In  this  case,  the  linear  relations  (3)  and  (4)  serves  as 
a  reasonable  starting  point  for  fine-tuning  or  optimization. 

2.4.3.  Implementation 

Some  preliminary  studies  of  incorporating  metabolic  information  into  the  IMRT 
inverse  planning  has  been  reported  by  our  group20  and  others  and  the  technical 
feasibility  of  planning  deliberately  non-uniform  dose  distributions  in  accordance  with 
functional  imaging  requirements  has  been  demonstrated.  In  our  preliminary  study,  a 
conventional  quadratic  objective  function  was  used  with  an  iterative  inverse  planning 
algorithm  for  the  optimization  of  the  system  with  inhomogeneous  dose  prescription 
specified  according  to  the  method  described  above.  The  generalized  quadratic 
objective  function  used  for  the  4D-dose  optimization  problem  reads  20 

F  =  *L  ir'Zrn  [Dc(n)-Dp(n)]2  ,  (5) 

where  Nc  represents  the  total  number  of  voxels  of  a  structure,  Dt(n)  is  the 

|alcuIaBd^dose“'S|oxel  n,  prescribed  dose  given  by.Eq.  (3)  or^4j 

dependii^^^yhS'if  the  voISJn;  belongs  to  the  target  or  normal  tissue^fThe 
weighting  factor  at  a  voxel  n  is  a  product  of  two  factors,  an  overall  factor  specific  to 
the  structure  o,  ra,  and  a  voxel  dependent  component  60,  rn,  describing  the  relative 
weighting  of  different  voxels  inside  the  structure. 
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For  testing  purpose  we  constructed  a  phantom  with  a  few  hypothetical 
metabolic  distributions  in  the  tumor  and  functional  importance  distributions  in  the 
sensitive  structure,  as  shown  in  top  row  of  figures  8.  The  use  of  phantom  case  with 
hypothetical  functional  data  allows  us  to  systematically  test  the  performance  of  the 
algorithm  effectively  without  going  into  the  technical  details  of  functional  imaging 
modalities.  In  figure  8  we  show  the  IMRT  plans  obtained  for  these  hypothetical 
situations.  Figure  9  shows  a  six-field  (0,  55,  135,  180,  225  and  305  degrees  in  IEC 
convention)  IMRT  glioma  case.  The  MRSI  metabolic  map  was  discretized  into  three 
discrete  levels  (figure  9A).  The  level  of  abnormality  at  a  point  is  characterized  by  an 
index  based  on  the  number  of  standard  deviations  (SD)  from  normal  values  of  the 
choline/NAA  ratio.  The  tumor  was  unifocal  and  44  Gy  was  prescribed  to  the  volume 
between  the  tumor  boundary  and  the  1st  abnormality  level  (AL=3)  and  64  Gy  was 
prescribed  to  the  highest  abnormal  region  (AL  between  5  and  7). 

2.4.4  Role  of  intra-structural  tradeoff 

Even  in  conventional  inverse  planning  scheme,  voxels  within  a  target  or  a 
sensitive  structure  volume  are  generally  not  equivalent  in  achieving  their  dosimetric 
goals  in  IMRT  planning.  Depending  on  the  patient’s  geometry,  beam  modality  and 
field  configuration,  some  regions  may  have  better  chance  to  meet  the  prescription 
than  others,  and  vise  versa.  It  has  shown  61  that  the  purposed  modulation  of  spatial 
penalty  distribution  is  more  advantageous  over  the  conventional  inverse  planning 
technique  with  structurally  uniform  importance  factors,  leading  to  significantly 
improved  IMRT  treatment  plans  that  would  otherwise  unattainable.  An  example  is 
given  in  figure  10,  in  which  the  isodose  curves  are  “pushed”  toward  the  target 
volume  and  the  dose  gradient  at  the  tumor  boundary  is  greatly  increased.  The 
significant  improvement  is  also  demonstrated  in  the  DVH  plots.  It  is  remarkable  that 
by  simply  modulating  the  spatial  importance  distribution  an  almost  uniform 
reduction  of  ~20%  (normalized  to  the  maximum  sensitive  structure  dose)  in  the 
sensitive  structure  dose  was  accomplished.  Conversely,  the  target  dose  can  often  be 
escalated  by  ~10%  while  keeping  the  radiation  toxicity  at  its  current  IMRT  level. 

The  intra-structural  tradeoff  plays  a  more  important  when  dealing  with 
biologically  heterogeneous  systems  since  non-uniform  dose  prescription  often 
aggregates  the  competitation  among  the  voxels.  The  approach  proposed  by  Shou  and 
Xing  can  be  easily  extended  for  the  determination  of  an  adequate  set  of  voxel 
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dependent  importance  factors  to  model  the  intra-structural  tradeoff.  Briefly,  once  the 
prescription  dose  is  given,  it  is  possible  to  quantify  the  degree  for  a  voxel  to  achieve 
its  dosimetric  goal  by  introducing  the  concept  of  dosimetric  capability  for  each  voxel 
in  a  target  or  sensitive  structure.  As  an  example,  in  figure  1 1  we  show  the  capability 
maps  obtained  for  a  uniform  dose  prescription  when  five  incicdent  beams  are 
invloved.  The  capability  of  a  voxel  represents  a  priori  dosimetric  knowledge  of  the 
system.  The  intra-organ  tradeoff  is  then  modulated  purposely  using  a  heuristic 
relationship  between  the  inherent  dosimetric  capability  and  the  voxel-based 
weighting  factors.  In  such  way,  we  can  impose  a  differential  penalty  scheme  and 
allows  the  system  to  suppress  potential  overdosing  spots  and  boost  the  potential 
underdosing  spots,  leading  to  a  solution  that  is  more  consistent  to  our  clinical 
expectation. 

The  voxel  based  penalty  scheme  can  also  be  used  as  a  means  to  fine-tune  the 
regional  doses.  Clinically,  it  happens  frequently  that,  after  optimization,  the  dose  at 
all  but  just  one  or  a  few  small  regions  are  satisfactory  and  thus  prevent  the  plan  from 
being  acceptable.  The  difficulty  is  that  the  location  of  the  hot/cold  region  in  inverse 
planning  is  generally  not  known  until  the  “optimal”  plan  is  obtained.  Consequently, 
an  “on-the-fly”  mechanism  is  highly  desirable  to  adaptively  fine-tune  the  dose 
distribution  after  a  solution  close  to  the  optimum  is  obtained.  Currently,  the 
modification  can  only  be  achieved  through  adjusting  structure  dependent  system 
parameters  (e.g.,  prescription,  importance  factors),  which  influence  not  only  the  dose 
at  the  region  of  interest  but  also  at  other  areas.  In  order  to  modify  the  dose  at  a 
specific  region,  in  principle,  one  can  use  ray-tracing  to  find  the  beamlets  that 
intercept  the  area  and  adjust  their  intensities  accordingly.  But  there  are  numerous 
ways  to  change  and  the  optimal  arrangement  of  the  beamlet  intensities  is  not  obvious. 
In  biologically  conformal  IMRT,  the  issue  becomes  more  urgent.  The  voxel 
dependent  penalty  scheme  provides  a  practical  solution  for  us  to  modify  the  local 
dosimetric  behavior  effectively,  as  has  been  demonstrated  in  recent  studies  of  our 
group  62,  63  Wu  et  al.  64 


2.4.5  Spectral  uncertainty 

In  practice,  molecular/functional  imaging  data  do  not  always  accurately  reflect  the 
actual  metabolic  level  over  the  entire  imaging  volume  because  of  some  technical 
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limitations  (in  MRSI,  for  example,  shimming  can  be  problematic  near  air-filled  cavities 
and  may  strongly  depend  on  the  surface  coil  SNR  on  the  spatial  position).  To  fully  utilize 
the  metabolic  information,  it  is  desirbale  to  develop  an  algorithm  to  numerically 
incorporate  the  spectral  uncertainties  (confidence  map)  into  IMRT  treatment  planning. 
A  statistical  analysis-based  inverse  planning  seems  to  be  ideally  suitable  for  this  purpose. 
Assuming  that  the  fluctuation  of  the  spectral  activity  or  the  prescribed  dose  Lf(n)  at 
voxel  n  is  specified  by  a  probability  distribution  Pn(Lf),  we  incorporate  the  P„(If)  by 
using  a  statistical  inference  technique.  Considering  that  currently  available  functional 
image  data  are  not  completely  reliable  and  that  missing  or  incomplete  spectral  data 
may  occurs  frequently,  such  type  of  technique  shuold  be  useful  to  minimize  the 
effect  and  generate  statistically  optimal  treatment  plans.  When  there  is  no  uncertainty 
in  the  spectral  data,  the  algorithm  reduces  to  the  conventional  inverse  planning  scheme. 


2.4.6  Biological  model  for  molecular/functional  image-guided  IMRT 

Dose-based,  and  more  recently,  clinical  knowledge-based  model  provides  an 
immediately  applicable  technique  for  generating  spatially  non-uniform  dose 
distributions.  However,  a  biological  model-based  approach  is  more  fundamental  and 
logical  in  dealing  biological  imaging  data  and  is  worth  of  a  detailed  investigation. 
Two  important  questions  in  the  biological  modeling  of  the  system  are;65  (i)  how  to 
determine  the  non-uniform  dose  prescription  provided  that  the  biology  distribution  is 
known;  and  (ii)  how  to  find  the  optimal  solution.  While  the  latter  problem  is  similar 
to  that  in  conventional  INRT  inverse  planning,  the  solution  to  the  first  problem 
entails  some  theoretical  considerations.  In  Sec.  2.4.2  we  have  used  the  metabolic 
abnormality  index  to  phenomenologically  characterize  the  tumor  burden.  Using 
radiobiological  model,  it  is  possible  to  relate  the  prescription  dose  to  the  more 
fundamental  radiobiology  parameters  to  optimize  the  cell  killing. 

Let  us  start  with  the  linear  quadratic  (LQ)  model.  We  include  the  effect  of  tumor 
cell  proliferation  but  ignore  the  quadratic  term.  The  model  parameters  include 
clonogen  density  (p),  radiosensitivety  (a),  and  proliferation  rate  (y).  The  time 
dependence  of  the  parameters  are  ignored.  The  tumor  control  probability,  TCPh  for  a 
tumor  voxel  i,  can  be  expressed  as 

TCP,  =  exp[-/?0,  V,  exp(-«, Di  +  y,  AT)] ,  (6) 


13 


where  V,  is  the  volume  of  voxel  i,  p0i ,  at  and  yt  represent  the  initial  clonogen 

density,  radiosensitivety  and  proliferation  rate  in  voxel  i,  respectively,  £>,  is  the  dose 
received  by  voxel  /,  and  AT  is  the  overall  treatment  time.  In  Eq.  (6),  yj  =  In  2  /  Tpl 

where  TPi  is  the  potential  doubling  time  in  voxel  i.  TCP  for  the  tumor  is  given  by 

TCP  =  Y[TCP, .  (7) 

/ 

A  constraint  from  the  normal  cells  within  the  tumor  volume  given  by 

ZmA=E,  (8) 

i 

should  be  applied  to  determine  the  tumor  dose  prescription,  where  m,  is  the  mass  of 
voxel  /,  E,  is  the  integral  dose  in  tumor.  The  problem  now  becomes  to  maximize  the 
TCP  under  the  constraint  of  equation  (6),  which  can  be  solved  using  the  method  of 

Lagrange  multipliers66.  When  the  mass  and  volume  are  equal  for  all  tumor  voxels, 
the  desired  prescription  dose  of  a  voxel  is  given  by 

D,  =^Dr~— (yr-ri)AT-±-  lnfe^l  (9) 

«/  «,  a,  atp0i  J 

IMRT  inverse  planning  optimization  can  proceed  by  numerically  maximize  the 
TCP  while  maintaining  the  NTCP  below  a  certain  limit.  One  can  also  take  a  “hybrid” 
approach  by  using  the  conventional  objective  function  with  the  above  dose 
prescription. 

2.4.7  Plan  review  tools 

The  sheer  volume  of  information  inherent  in  3D  treatment  designs  and  the 
corresponding  dose  distributions  make  display  and  objective  assessment  problematic. 
Details  of  a  dose-distribution's  spatial  characteristics  can  be  obtained  by  examining  2D 
isodose  curves  in  a  slice-by-slice  fashion;  however,  this  is  a  quasi-quantitative,  time- 
consuming  process  and  is  not  an  efficient  way  to  compare  competing  plans  even  for 
conventional  IMRT.  In  the  presence  of  an  additional  degree  of  freedom  (metabolic 
abnormality),  the  problem  is  exacerbated  by  the  breakdown  of  uniform  dose  assumption 
within  the  target  volume.  One  of  the  commonly  used  approaches  is  the  reliance  on  data 
reduction  techniques  in  the  quantitative  assessment  of  alternative  plans.  DVH  is  one  of 
the  most  widely  used  data  reduction  techniques.  This  technique  enables  the  ready 
reduction  of  the  complex  3D  data  set  of  a  treatment  design  into  the  2D  display  of  the 

fractional  volume  of  a  given  structure  receiving  doses  within  a  particular  range. 
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Unfortunately,  this  tool  becomes  invalid  for  4D  plan  evaluation  because  of  possibly  non- 
uniform  biological  status  of  the  involved  structures.  4D  IMRT  technique  requires  new 
plan  review  tools  to  facilitate  the  quantitative  comparison  of  plans.  The  following  are  a 
few  tools  that  are  potentially  useful  for  the  plan  review  of  molecular/functional  image 
guided  IMRT  plans. 

a) .  Effective  isodose  dose  distribution:  The  effective  dose  at  a  voxel  is  defined  as  the 
ratio  of  the  physical  dose  and  the  prescribed  dose.  This  distribution  considers  both  the 
spatial  dose  distribution  and  the  metabolic  map  and  provides  intuitive  information  on  the 
geometric  location  of  underdosing  or  overdosing  regions.  In  this  way,  we  can  use 
conventional  wisdom  to  evaluate  a  4D-dose  distribution.  The  DVH  corresponding  to  the 
effective  dose  distribution  is  also  useful. 

b) .  DVH  clusters:  In  practice,  not  all  underdosing/overdosing  are  equally  significant 
and  underdosing/overdosing  at  a  certain  metabolic  level  maybe  more  acceptable  than  at 
other  metabolic  levels.  A  cluster  of  DVHs,  each  corresponding  to  an  incremental  range 
of  metabolic  activity  of  interest,  may  provide  useful  tool  to  address  the  issue.  The  cluster 
of  DVHs  can  be  used  to  check  the  overall  dosimetric  behavior  at  an  individual  metabolic 
level.  Figure  1 1  represents  an  example  of  a  3-level  DVH  cluster.  For  a  sensitive  structure 
with  functional  data  available,  similar  technique  applies. 

c) .  Functional  dose-volume  histogram  (FDVH):  Distribution  of  functional  importance 
appears  to  be  heterogeneous  in  some  normal  organs  and  functional  imaging  modalities 
such  as  MRSI  or  PET/SPECT  may  provide  valuable  information  about  the  spatial 
distribution  of  the  functional  importance.  The  FDVH,  originally  introduced  by  Lu  et  al. 

67,  Marks  et  al.  68,  and  Alber  and  Nusslin  69,  may  prove  to  be  a  useful  plan  review  tools. 
A  similar  histogram  function  can  be  introduced  for  the  tumor,  but  its  usefulness  needs  to 
be  justified. 

d) .  Modified  TCP  and  NTCP  calculation  tools:  The  conventional  TCP  and  NTCP 
formula  67  50,  70, 71  nee(j  t0  be  modified  to  take  into  account  the  heterogeneous  biology 
distribution  72-74  This  modification  should  be  straightforward  if  the  spatial  distributions 
of  radiobiological  parameters  are  known.  Although  it  is  difficult  to  obtain  quantitative 
results  from  the  model  calculation  because  of  the  uncertainties  in  the  parameters, 
qualitative  conclusions  regarding  the  deliberately  non-uniform  irradiation  scheme  can  be 
drawn  and  may  shed  useful  insight  into  the  problem  74, 


3.  Conclusion 

The  success  of  radiotherapy  critically  depends  on  the  imaging  modality  used  for 
treatment  planning  and  the  level  of  integration  of  the  available  imaging  information. 
The  use  of  functional/metabolic  imaging  provides  us  much  more  than  a  tool  to  better 
delineate  the  boundary  of  a  tumor  target.  Together  with  anatomical  CT  or  MRI 
images,  functional  imaging  affords  valuable  4D  data  (3D  structural  plus  ID 
metabolic)  for  both  tumor  and  sensitive  structures,  valuable  for  guiding  us  to  design 
spatially  non-uniform  dose  distributions  to  deliver  high  doses  to  where  the  tumor 
burdens  are  high  and  differentially  spare  the  sensitive  structures  according  to  the 
functional  importance  distributions.  The  integration  and  utilization  of  the  functional 
data  in  radiation  therapy  treatment  planning  become  increasingly  important  to 
improve  clinical  cancer  management.  While  it  is  straightforward  to  modify  the 
radiation  portals  to  accommodate  any  changes  in  treatment  volume,  new  methods  of 
dose  optimization  and  medical  decision-making  must  be  developed  to  take  full 
advantage  of  the  metabolic  information  and  IMRT.  How  to  achieve  biologically 
conformal  doses,  instead  of  the  geometrically  conformal  dose  distribution,  presents  a 
new  challenge  to  radiation  oncology  discipline.  Hopefully,  with  the  efforts  from 
muti-institutions,  the  new  approach  of  imaging,  planning  and  decision-making  will 
be  resolved.  Ultimately,  whether  using  deliberately  inhomogeneous  dose  distributions 
obtained  under  the  guidance  of  functional  imaging  such  as  MRSI  can  improve  patient 
survival  and  reduce  the  side  effects  associated  with  radiation  treatment  should  be 
established  through  extensive  clinical  trials. 
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